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SUMMARY 


A  numerical  discrete -element  method  of  wave  motion  analysis  is 
summarized  and  extended  for  problems  Involving  Infinite  or  semi -Inf Ini te 
solid  media  in  plane  and  axl-symmetric  conditions.  Space  discretization  of 
i  solid  medium  is  accomplished  through  a  lumped-parameter  discrete-element 
model  of  the  medium,  whereas  the  time  discretization  is  embedded  within  a 
general  numerical  integrator.  This  invariably  leads  to  a  system  of  finite 
difference  equations;  thus,  the  required  mathematical  conditions  for  numerical 
stability  can  be  developed  on  the  basis  of  available  finite  difference  theory. 
Explicit  stability  conditions  for  plane  and  axi-symmetric  problems  are  presented. 

Calculations  of  wave  motions  in  an  Infinite  or  semi-infinite  space 
can  be  confined  to  a  finite  region  or  Interest  if  the  region  is  terminated  by 
suitable  "transmitting  boundaries"  such  that  no  significant  reflections  are 
generated  at  these  artificial  boundaries.  Based  on  the  concept  of  a  step-wise 
transmission  of  D'Alembert  forces,  a  general  transmitting  boundary  was  developed 
for  the  discrete-element  method  of  analysis.  The  boundary  was  verified 
extensively  through  actual  calculations  of  plane  strain  and  axl-symmetric 
problems,  including  those  with  layered  half-spaces,  elastic-plastic  systems, 
and  a  problem  involving  long  calculation  time. 
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I.  INTRODUCTION 


I  .  1  Statement  of  Problem 

The  prediction  of  wave  motions  in  solid  media  often  requires  discrete 
numerical  solution  techniques.  One  such  method  is  the  calculus  of  finite 
differences.  The  usefulness  of  the  method  has  been  enhanced  greatly  through  the 
development  of  discrete-element  lumped-parameter  models  of  continuous  media, 
through  these  discrete  models  the  required  difference  equations  (usually  centered 
difference)  can  be  derived  directly  on  the  basis  of  clear  physical  considerations. 
This  removes  many  of  the  difficulties  usually  associated  with  boundary  conditions 
when  the  necessary  discretizations  are  imposed  or  applied  purely  through  mathematical 
operations  as  is  ordinarily  dene  with  the  application  of  finite  difference  equations. 

In  any  discrete  method  of  analysis  (including  finite  element  [13] )»  the 
space  domain  must  necessarily  be  limited  to  some  finite  region.  When  such  discrete 
methods  are  applied  to  wave  propagation  problems  involving  extensive  space  domains 
(could  be  infinite  or  semi  finite),  the  amount  of  calculations  required  may 
conceivably  become  excessive,  to  the  extent  that  it  becomes  uneconomical  ,  or  even 
impossible,  to  perform  the  calculations.  This  limitation  of  the  discrete  formulation 
may  be  removed  if  the  required  calculations  can  be  confined  to  a  finite  space  domain 
of  interest,  such  that  there  are  no  reflections  from  the  terminating  boundaries  of 
the  finite  region.  Unfortunately,  there  are  no  such  boundaries  that  are  generally 
available  for  this  purpose. 

The  objective  of  this  study  is  to  develop  a  "transmitting  boundary"  that 
can  be  used  for  numerical  wave  propagation  calculations  in  plane  and  ax i -symme tr i c 
med i a . 

I . 2  General  Solution  Method 

The  transmitting  boundary  developed  herein  is  intended  to  be  used  with 
a  general  method  of  numerical  solution  that  is  essentially  based  on  the  mathematical 
theory  of  finite  differences.  The  technique  of  formulation  requires  discretization 
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of  the  space  domain  as  well  as  the  time  domain.  Proper  bases  for  such  discreti- 
zntions  can  be  accomplished  through  a  physical  discretization  of  the  space  domain 
through  the  use  of  discrete  element  models.  It  should  be  emphasized,  however, 
that  for  dynamic  problems  of  wave  propagation,  a  formulation  of  the  problem  based 
merely  on  a  proper  discretization  of  the  space  domain  does  not  necessarily  lead  to 
a  solution  (even  in  an  approximate  sense).  To  assure  a  valid  approximate  solution, 
certain  mathematical  requirements  must  be  satisfied. 

In  the  method  of  formulation  suggested  herein,  the  discretization  of  the 
space  domain  is  accomplished  through  the  use  of  appropriate  mathematically 
consistent  d i scrate-e I  ament  models  of  the  solid  media,  whereas  the  time  discreti¬ 
zation  is  embedded  within  the  numerical  integrator.  The  space-time  discretization 
is  mathematically  equivalent  to  a  finite  difference  formulation;  consequently,  the 
mathematical  relations  required  to  assure  stability  and  convergence  of  the  resulting 
solutions  can  be  studied  from  the  standpoint  of  the  stability  of  finite  differences. 
Explicit  stability  conditions  were  obtained  for  the  algorithmic  schemes  suggested 
he  re  i  n . 

1 .3  Scope  of  Investigation 

The  general  method  of  analysis  for  wave  propagation  has  previously  been 
applied  extensively  to  numerous  problems  of  solid  media,  including  large-scale 
problems  of  wave  and  ground  shock  predictions  from  blast  loadings,  and  earthquake- 
induced  ground  motions.  A  comprehensive  discussion  of  the  method  of  formulation 
and  solution  is  summarized  in  Chapter  II,  including  new  mathemat i cal  results 
pertaining  to  cal  cu!  ational  stability  requirements. 

The  main  scope  of  this  study  pertains  to  the  development  of  a  transmitting 
boundary.  The  major  concepts  involved  in  the  various  developments  are  presented  in 
Chapter  III.  Certain  aspects  of  the  developments  require  semi-empirical  investigations 
as  explained  in  Chapter  III;  the  concepts  presented  herein  represent  the  final 
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results  of  numerous  schemes  that  were  considered  and  tested  numerically  in  these 
investigative  studies,  some  of  these  are  illustrated  in  Chapter  IV. 

The  results  of  the  specific  investigations  may  be  summarized  as  follows 

(i)  For  the  one-dimensional  plane  elastic  propagation,  the  step-wise 
transmission  of  D'Alembert  forces  is  theoretically  exact;  this  is  verified  by  the 
specific  numerical  calculations. 

(ii)  For  the  two-dimensional  propagation  under  plane  strain  and  axi- 
symmetric  conditions,  the  exact  speed  of  transmission  of  the  D'Alembert  forces 
is  not  known;  certain  combination  of  the  dilatational  and  shear  velocities  of  the 
material  was  found  to  be  suitable.  For  low  velocity  material,  this  was  found  to 
give  excellent  results;  however,  for  fast  material  (i.e.,  c^  <  3>000  fps)  certain 
high-frequency  oscillations  were  observed  on  unloading.  Such  oscillations, 
however,  can  be  corrected  with  the  use  of  artificial  viscosity,  which  is  normally 
required  also  for  numerical  stability  in  axi-symmetrlc  calculations.  The  trans¬ 
mitting  boundary  is  applicable  to  problems  with  layerings,  as  well  as  general 
loadings  including  periodic  loading  histories. 

(iii)  Extension  of  the  transmitting  boundary  to  elastic-plastic  media 
requires  the  inclusion  of  the  fully  plastic  velocity  of  propagation,  as  well  as 
the  propagation  of  the  elastic  precursor  waves. 

A  number  of  schemes  using  different  combinations  of  elastic  and  plastic 
velocities  were  devised  and  tested;  these  led  to  the  development  of  an  elastic- 
plastic  scheme  for  the  transmitting  boundary.  Although  not  fully  verified  for 
all  cases,  the  algorithm  presented  herein  is  suitable  for  the  elast ic-perfect ly 
plastic  problems  examined  herein,  including  those  in  which  plastic  waves  are 
involved  at  the  terminating  boundary. 

The  Investigation  included  also  an  exploratory  application  of  the 
transmitting  boundary  to  a  finite-element  type  of  formulation.  The  results. 
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although  not  very  conclusive,  indicated  that  the  step-wise  transmission  of 
D'Alembert  forces,  which  is  the  basis  of  the  proposed  transmitting  boundary, 
may  not  be  directly  applicable  to  the  finite  element  method  of  analysis. 
Modifications  of  the  basic  concept,  however,  may  make  it  applicable;  additional 
studies  for  this  purpose,  however,  are  required. 

Finally,  an  alternate  finite  difference  scheme  was  developed  to  serve 
the  same  purpose  as  that  of  the  transmitting  boundary.  This  is  based  on  certain 
mathematical  properties  of  wave  propagation,  and  can  be  used  in  place  of  an 
explicit  transmitting  boundary;  however,  the  scheme  is  limited  so  far  to  purely 
elastic  calculations. 

I  .4  Ac  know  I  edeynen  ts 

The  work  ruportad  herein  was  conducted  for  the  Defense  Atomic  Support 
Agency  as  Task  I  of  Contract  D ASA-0 I -69-0040. 

A  number  of  other  technical  personnel  were  involved  In  conducting  the 
research  studies.  Specifically,  Or.  L.  A.  Lopez  and  Dr.  A.  H.  M.  Sameh  assisted 
in  the  development  of  the  requisite  computer  programs,  and  Dr.  Y.  Uckan  assisted 
in  the  development  of  some  of  the  mathematical  relationships. 

I . 5  Notations 

The  following  list  of  symbols  and  notations  is  used  in  this  report: 
velocity  of  dilatational  waves  in  an  elastic  medium 
cg  velocity  of  shear  waves  in  an  elastic  medium 

c  wave  velocity,  in  general 

.  I 

E  vector  representing  solution  error 

I  direct  strain 

G  Lame  constant  (shear  modulus) 

h  unit  dilatational  transit  time  h  »  Ay/c. 

0 

k  unit  shear  transit  time  k  «  Ay/cs 
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•  >* 


K 

*vS 

I 

i 

tmdginary  constant ,  '/-I 

Index  for  r  (or  x)  coordinate  of  a  point 

j 

index  for  z  (or  y)  coordinate  of  a  point 

n 

index  for  time  level  *  t/  At 

p 

index  for  distance  of  a  mass  point  from  the  r-  (or  x-)  axis 

-  r/  Ar 

q 

index  for  distance  of  a  mass  point  from  the  z-  (or  y-)  axis 

•  z/  Az 

V  qz 

artificial  pseudo-vi scous  stress  In  r-  and  z-di rections, 

respectively 

r.  •,  i 

cylindrical  coordinate  system 

rl 

distance  of  the  i  th  point  from  axis  of  symmetry 

t 

time  variable 

u 

displacement  of  a  mass  point  In  r-  (or  x-)  direction 

u 

velocity  of  a  mass  point  In  r-  (or  x-)  direction 

u 

acceleration  of  a  mass  point  in  r-  (or  x-)  direction 

u 

dependent  variable  vector 

V 

Fourier  coefficient  vector 

w 

displacement  of  a  mass  point  In  z-  (or  y-)  direction 

X 

£ 

V 

an  indicator  of  plastic  yielding 

*»  y 

y 

Cartesian  coordinate  system  in  a  plane 

non-Cartesian  local  reference  system  in  a  plane 

or 

angle  between  x-  and  x-axes;  also  as  given  in  Eq.  (3.24) 

e 

a  parameter  in  Newmark's  B*  method;  also  as  given  in  Eq.  (3.25) 

r 

coefficient  of  li.-tear  artificial  viscosity 

5 


V 


av 


Ar,  ae,  a z 
At 

Ax,  Ay 

X 

V 

t 

mi 

0 

a 


a’ 


shearing  strain 

volume  of  an  element  of  a  mass  point  of  the  lumped-parameter  model 
uniform  mesh  sizes  In  r-,  6-,  and  z-di recti ons,  respectively 
time  increment 

uniform  mesh  sizes  in  x-  and  y-di rections,  respectively 
Lame  constant  of  medium 
Poisson  ratio  of  medium 

factor  representing  error  growth  from  one  time  step  to  the  next 

mass  density  of  medium 

direct  stress 

shearing  stress 

angle  between  y-  and  x-axes 

ratio  of  space  meshes,  Q  *  Lr/  Az 
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.  DISCRETE  FORMULAT 'ON  OF  SPACE-TIME  PROBLEMS 

Analytic  formulation  of  space-time  problems  of  solid  media  invariably 
leads  to  a  system  of  partial  differential  or  in tegro-d i f fe ren t I al  equations. 

The  solutions  of  these  equations  are  the  relevant  physical  quantities  desired  in 
a  problem,  which  are  generally  continuous  functions  of  space  and  time.  Such 
continuous  analytic  solutions,  however,  are  quite  difficult  to  find,  except  for 
certain  simple  problems.  For  problems  of  practical  engineering  significance, 
continuous  solutions  are  often  virtually  impossible  to  obtain;  in  these  instances, 
approximate  numerical  solutions  are  the  only  practical  alternative.  Invariably, 
such  numerical  solutions  represent  discrete  approximations  to  the  continuous 
functions  of  space  and  time.  This  is  often  accomplished  by  formulating  an 
associated  set  of  difference  equations,  the  corresponding  solutions  of  which 
yield  a  discrete  function  of  space  and  time  --  i .e. ,  it  is  a  function  defined 
only  at  discrete  points  in  space  and  at  discrete  instants  of  time. 

For  space-time  problems,  a  discratized  formulation  does  not  necessarily 
imply  a  solution;  i.e.,  the  solution  of  the  discretized  system  of  equations  may  or 
may  not  lead  to  a  valid  approximation  of  the  "correct"  solution.  Certain  mathema¬ 
tical  requirements  must  be  satisfied  in  the  numerical  solution  process  to  assure 
the  validity  of  the  approximate  solution.  These  include  questions  of  convergence 
and  stability  of  a  particular  numerical  scheme. 

2. 1  Space  Discretization 

A  discrete  formulation  of  a  space-time  problem  requires  discretization 
of  the  space  domain,  as  well  as  the  time  domain.  It  will  be  seen  later  that  the 
resulting  solution  will  depend  on  the  relation  between  the  space  and  time 
d i sere t i zations . 

Discretization  of  the  space  domain  may  be  done  mathematically  in  terms 
of  finite  differences;  alternatively,  the  space  discretization  may  be  accomplished 


7 


on  a  physical  basis.  This  latter  approach  requires  I umped-parameter  or  discrete* 
element  (e.g.,  finite  element)  models  of  the  otherwise  continuous  space  problem. 

In  any  such  models,  the  basic  physical  quantities  are  defined  appropriately  at  a 
finite  number  of  discrete  positions  in  the  space  domain.  The  resulting  set  of 
equations  formulated  for  the  discrete-element  model  may  also  be  (centered)  finite 
difference  equations;  such  a  property  may  be  referred  to  as  "mathematical  consistency" 
of  the  discrete-element  model. 

Lumped-parameter  models  for  two  and  three  dimensional  solid  media  have 
been  developed  and  described  extensively  elsewhere  [  I,  2,  4).  However,  for  the 
sake  of  coherence  and  completeness,  the  models  for  plane  strain  and  axi-symmetr.c 
media  are  summarized  herein. 

Plane  Strain  Model  - —  In  rectangular  Cartesian  coordinates,  the  plane 
strain  model  is  shown  graphically  in  Fig. 2.1.  As  are  cannon  with  models  of  the 
same  type,  the  basic  elements  of  the  plane  strain  model  consist  simply  of  mass 
points  and  stress  points,  at  which  the  particle  motion'  (accelerations,  velocities, 
and  displacements)  and  stress-related  quantities  (such  as  stresses  and  strains)  of 
the  solid  medium  are  defined,  respectively. 

Based  on  the  fundamental  principles  of  dynamics,  the  equations  describing 
the  particle  motions  of  the  mass  points  in  the  model  are  (given  for  mass  point 
( i  ,  j )  at  t ime  t) : 


•i<i+u) - ®i(M.j>  ** (i.j+0 - t‘o.j-i>  . 

+  _*2 - - BUt(l  ,j  ) 


Ax 


Ay 


and  • 


<r!(i.j+D  -  •Jo.j-i)  (ih.j)  -  r* v(i-i,j) 

J. - +  -JU£ - - - S* - pvt(I  ,j) 


Ay 


Ax 


(2.1) 


(2.2) 


which  are  clearly,  respectively,  the  usual  centered  f ini te  dl fference  analogues 
of  the  following  differential  equations  for  plane  strain  motions: 
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(2.3) 


!!x , 

dTxv 

a2u 

ay 

— £  .  —21 

a! v 

=*  0  2 

*y 

*x 

at2 

Eqs.  (2.1)  and  (2.2)  are  quite  general  and  apply  to  any  material  in  plane 
strain,  Including  materials  with  non  I i near- ine last ic  behavior.  However,  for 
problems  where  small  strains  are  Involved,  the  nontrivial  strain-displacement 
relations  for  plane  strain  condition  are: 


*x  (i *J  ) 


u(i‘l.i)  ~  u(i-l.i) 

Ax 


«y  »J  > 


Vxy(t,j) 


v  ( i ,  j  +1 )  -  v (I  , j - 1 ) 
tv 

u(i,jH)  -  u(l,J-l)  v(l*l,j)  -  v  (I- 1 ,  j  ) 

- - - -  ♦  -  - - 

tv  Ax 


(2.5) 


and  the  corresponding  strain-rates  can  be  obtained  by  simply  replacing  the 
displacements  In  Eq . (2.5)  by  the  corresponding  velocities. 

Through  the  stress-strain  relationships  (or  constitutive  equations) 
of  the  material,  it  is  clear  that  Eqs. (2.1)  and  (2.2}  can  be  expressed  in  terms  of 
the  particle  motions;  specifically,  in  terms  of  the  displacements  u,  v  and 
velocities  u,  v  (or  Au,  Av).  In  particular,  if  the  material  is  linearly  Hookean, 
then  Eqs.  (2.1)  and  (2.2),  respectively,  become 
2 

u  ( I  —2 ,  J  )  -  2u(l,J)  +  u  (1  +2  ,J )  I 

(Ax)2  L 

c2 

+  — ^  [  u  (I  ,j-2)  -  2u  (i ,  j  )  +  u  ( i , j  +2 )  ] 

(Ay)2  L 
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2  _  2 

+  — - -  ["  v(i+l  ,j+l)  -  v(i+l ,j-l)  -  v  (  i  — 1  ,j+l)  +  v(i-l.j-l)]  =  u(i  ,  j  )  (2.b) 

Ax*  Ay 


and  , 

2 

— I  v(i,j-2)  -  2v(i,.i)  +  v  ( i  ,  i+2)  1 
(Ay)2  L 

c2 

+  — v ( i  -2  , j  )  -  2v(  i  ,j  )  +  v(  i+2  ,j  )  1 
(Ax)*  L 


2 

Cd  " 


Ax*  Ay 


I"  u(i+l  ,j  +  l )  -  u(  i+l  ,j-l )  -  u(i-1  ,j+l )  +  u ( I -1  , 1-1 )  J  *  v(i  ,j  )  (2.7) 


Plane  Strain  Model  in  Non -Cartes i an  Reference  ---  I n  model  I ng  sol  I d 
media  with  irregular  boundaries,  or  containing  irregular  openings  and  Inclusions, 
the  discretization  of  the  space  domain  must  include  such  geometric  irregularities. 
For  these  purposes,  the  discretization  of  the  space  domain  may  be  described  in 
non-Cartesian  (i.e.,  nonorthoyonal )  coordinates.  Fig.  2.2  shows  such  a  description 
for  the  plane  model.  In  this  general  frame  of  reference,  a  local  non-Cartesian 
coordinate  and  a  global  Cartesian  coordinate  systems  are  required.  The  particle 
motions  of  the  mass  points  are  referred  to  the  global  system;  however,  at  a  stress 
point  which  is  incident  on  four  neighboring  mass  points  (see  Fig.  2.2),  the  strains 
are  intrinsically  defined  in  the  local  reference.  A  transformation  between  the 
local  references  (generally  different  from  stress  point  to  stress  point)  and  the 
global  Cartesian  reference,  therefore,  is  required. 

Denoting  the  global  coordinates  as  x,  y  and  the  local  coordinates  as 
x,  y,  the  required  transformation  for  the  plane  model  Is  given  by, 

x  -  -  d,)(x  sin  in  -  y  cos  «)  (2.8a) 

y  „  __J -  (x  sin  Q  +  y  cos  a)  (2,Rb) 

sin(<*  -  tr) 
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where  a  and  jj  are  the  angles  between  the  x-axIs  and  the  x-  and  y-axes , 
res pec  Li vel y ,  as  shown  in  Fig.  2.2a. 


The  strain  components  at  a  general  stress  point  (assuming  small  strain 


theory) 

then 

become , 

t 

8u 

8u 

8  x 

X 

8x 

8x 

8x 

8y  8y 

8v 

8  v 

8x 

+  *2.  .  8£ 

*y 

8  y 

8x 

8  x 

8y  8y 

V 

_ 

+  5* 

8u 

8x  |  8u  8y 

5  v 

■  • 

8  x 

8y 

8x 

~  a?. 

8y  8  y  8y 

8  x 

8x 

i n  wh i ch 

the 

parti 

i  al  der i vat i ves 

8x 

Ax 

’  8y* 

a£ 

8x 

(2.9a) 

5  v  8  y 

8y  8x 

and  ^  are  obtained  from  Eq .  (2.8) 


and  for  stress  point  0  in  Fig.  2.2a  which  Is  Incident  on  mass  points  I,  2,  3  and  k 


r  *0  = 

Jl! 

*  "»  ;  ( 

U2  - 

^  8x  -0 

Ax 

8y  '0 

Ay 

r  av  = 

-  v,  . 

->  .  ( 

'_8v-\  = 

v2  ~ 

_^4 

^  8x  4) 

Ax 

-  J 

'  8y  y0 

Ay 

where 

Ax  and  Ay 

are  the  mesh 

1 engths 

in  the 

local  reference  system.  Using 

these 

re! at i ons  , 

the 

strains  at 

stress 

point  0 

,  therefore,  are: 

e  = 


u2  '  u4 


1  r  Ul  '  U3  . 

- (  - - —  sin  uj  -  - — 

sin  {'u-a)  Ax  Ay 


s  i  n 


o 


1 


E  = 

y 


s  i  n  ( urar)  Ax 


( 


V1  '  V3  v2  v4 

cos  m  + 


cos 


Ay 


O 


xy 


c 


sin(urw)  Ax 


U1  ’  u3  u2  u4 

COS  (JU  +  -  COS  (Y 


Ay 


V1  •  V3 


Ax 


s*n  u) 


v2  '  v4 

Ay 


s  i  n 


°o 


(2.9b) 
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for  a  specified  stress-strain  equation,  the  associated  strsss  cornpon.  n's 
at  stress  point  0  can  then  be  detcrnined  wi  th  the  strains  given  above;  thes  • 
stresses  are  defined  in  the  global  (x,y)  reference.  By  suitable  tensorial 
irans format  ion  [7,  11  ],  the  stresses  in  the  local  reference  can  bj  obtained, 
which  are: 


I  2  2 

a-  =  1  ■■■  -  (0  sin  it  *  t  sln2tu+  0  cos  m  ) 

X  sin(«nv)  x  xy  y 


I  ?  2 

a-  =  — —  —  ■■  (o  sin'o  -  t  sin  2or  +  o  cos  a  ) 

y  .  t  \  .<  xy  y 

sinu'-o) 


(2.10) 


T._  = 


sin(,p-ar) 


j  -o  sh  a  sin  w  +  t  sin(t*f»)  -  ®  cos  or  cos  » 

L  ^  xy  y  -j 


The  stresses  at  all  the  stress  points  surrounding  a  mass  point  must  be  determined; 
having  these  stresses,  the  corresponding  forces  acting  on  the  mass  point  can  be 
determined  in  the  appropriate  local  reference.  These  forces  may  be  transformed 
into  the  global  refarence,  thus  permitting  the  writing  of  the  equations  of  motions 
of  the  mass  point  in  question,  yielding  for  mass  point  0  (see  Fig.  2.2b), 


I  -x.  ""  *1  “  V,  =OS  Wl)  *  ^3  (%  "  Txy3  co‘ 


Ay,  (•  si 
xl 


+  A*2  <-ax2  s,"»  ®2  +  Txy2  COS  °2)  +  A*4  {\  *4  *  \y4  **•  *4> 


*P  AV  Uq 


(2.11) 


and , 


Ax2  (Cy2  C0S  a2  *  Txy2  ®2>  +  A*4  (\y4  *ln  *4  *  \  CO»  V 
+  A*1  (t  sin  *,  -  •  cos  «,)  +  A? 3  (a  cos  *3  -  r  sin  *3) 
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2«>  AV  v0 


(2.12) 


where  AVq  is  the  elemental  volume  [  7,  II  }  of  nets  point  0. 

Using  the  straln-displacenent  equations  of  Eq.  (2.9b)  and  Hooke's 
stress-strain  equations,  the  equations  of  motions,  Eqs.  (2.11)  and  (2.12),  become 
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and , 
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u8 
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(2.14) 

for  *  -  0  and  i  •  «/2  (orthogonal  coordinates),  Eqs.  (2.13)  and  (2.14)  reduce  to 
Eqs.  (2.6)  and  (2.7),  respectively. 

Where : 

^  -  Ay  /  AS 

(X  ♦  c)  sin2*  ♦  6  (X  ♦  G)  sin  *  cos  a 

sin  (x  •  p)  sin  (#  •  j  ) 
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c  = 


l  * 


H  = 


J  «= 


(  X*  G)  sfn  x  sin  y  *  G  cos  (uj  -  or) 


sin  (ul-  -  nr) 


D 


X  sin  to  cos  or  ♦  G  sin  ry  cos  w 
sin  (tc  -  or) 


X  sin  »  cos  (j)  +  G  cos  a 


rin  (uj  -  a) 


F  - 


sin  (oj  -  or) 


2 

(X  ♦  G)  sin  y  cos  n  (X  *■  G)  cos  ud  +  G 

sin  (x>  -  or)  sfn  -  a) 

2 

(X  +  G)  cos  or  cos  gj  +  G  cos  (tu  -  <*)  (X  +  G)  cos  a  +  G 

-  ;  «  - - 

sin  (iu  -  or)  sin  (<u  -  or) 


Axl-symmetrlc  Model  ---  In  cylindrical  Cartesian  coordinates,  the 
axl -symmetric  model  shown  In  Fig. 2.3  can  be  used  to  formulate  problems  of  axl- 
symmetrlc  wave  motions.  Generally,  In  this  case,  two  sections  of  the  model  are 
required  to  describe  the  model,  and  the  corresponding  equations  of  motions  In 
terms  of  stresses  are. 


q*(l+i,J)  -  cr*(l-i,J)  t*2(i,J+0  -  t*20,J-0 


Az 


Az 


+  - * -  -  p  U  (I, 

r  (I  ,J) 


J) 


and. 


qjO  J+0  -  q2(l,J-l)  +  t ( l l , J )  -  t*z(I-I,J) 


Az 


Ar 


T;2d,j) 


r(l,J) 


o  (I,J) 


(2.15) 


(2.16) 


In  this  case,  the  strain-displacement  relations  are: 
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u(l+l,j)  -  U  ( I  -  I  ,  j  ) 


« 
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e 

z 


c 


e 


Y 


rf) 


Ar 

w(t  ,j  +  l )  -  w(i ,  j  - 1 ) 

Az 

u  (• ,] ) 
r  ( I ,  j  ) 

u  ( i  ,  j  + 1 )  -  u(t,j-l)  w(l+1,J)  -  w  (1  - 1 ,  j  ) 

Ai  +  Ar 


(2.17) 


Again,  if  the  material  equations  are  specified,  then  through  these  strain- 
displacement  relations  (plus  the  corresponding  strain  rate-velocity  relations), 
Eqs .  (2.15 )  and  (2.16)  can  be  expressed  In  terms  of  u,  v  and  u,  v  (or  Au ,  Av).  For  a 
Hookean  material,  Eqs.  (2.15)  and  (2.16)  become 


X  +  2G  p  ■.  X  +  2G  p  -i 

[  u(i-2,J)  -  2u  (t ,  j )  +  u(l+2,J)  J  +  -  ~  ■  -  [  u(i+l,J)  -  u(l-),j)  J 


(Ar ) 


r(Ar) 


X  +  2G 


™u(I,j)  +  — 5~r  u  (I ,  j-2)  -  2u  (I  ,j )  +  u  (l ,  j  +2)  1 
2  (Az)2  1  J 


X  *•  G 


Ar  •  Az 


[w  (I  +1  ,j-»l)  -  w(?+l,J-l)  -  w(l-l,j+l)  +  w(l-l,j-l)]  =  p  U(!,j)  (2.18a) 


and , 


X  +  2G 


(Az)' 


[w(l  ,j-2)  -  2w  (l ,  j  )  +  w(i  ,J+2)  J  +  —  [w(!+l  ,J  )  -  w(?-l  ,J  )  J 


X  +  G 


+  - T  j*w(i-2,j)  -  2w(i,j)  +w(l+2.j)l  +  -  [udJ+D  -  u(l,j-1)l 

(Ar)2  L  J  r-Az  L  J 


X  +  G 


Ar-Az 


[u(l+l,J+1)  -  u(l-l,J+l)  -  u(l+l,j-l)  +  u  (1-1  ,J-1 )]  *  p  w(!  ,j )  (2.18b) 
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2  .2  Time  Discretization 


In  each  of  the  cases  described  above,  we  have  a  set  of  differential- 
difference  equations;  e.g.,  Eqs.  (2.1)  and  (2.2)  for  the  plane  strain  case.  We 
observe  that,  so  far,  the  discretization  has  been  confined  only  to  the  space 
domain;  no  discretization  of  the  time  domain  has  been  Introduced,  and  presumably 
time  may  be  continuous.  However,  in  order  to  obtain  solutions  to  problems  involving 
general  time  functions,  discretization  of  the  time  domain  wilt  normally  also  be 
necessary.  The  resulting  solutions  will,  accordingly,  be  defined  only  at  the 
various  discrete  time  instants.  Such  solutions  are  usually  obtained  through  a 
step-by-step  numerical  integration  procedure.  There  are  numerous  integration 
schemes  available  for  this  purpose.  However,  for  a  large  class  of  second-order 
systems  encountered  In  wave  propagation  problems,  a  general  and  convenient  method  for 
this  purpose  is  the  Newmark  0-lntegrator  [9].  This  Is  a  step-wise  recursive  method 
of  numerical  quadrature.  The  basic  recursive  relations  for  advancing  a  small  but 
finite  time  step  At  are  the  following: 


u(t+At)  -  u(t)  +  iAt  t  u  (t )  +  U  (t+At)  ] 
v  (t+At )  =  v  (t )  +  iAt  [  v  (t )  +  v(t+At)  ]  > 

u  (t+At)  -  u  (t )  +  At-u(t)  +  (At)2  (i-B)  u(t)  +  (At  )2p  u(t+At) 

v (t+At)  =>  v(t)  +  At’v'(t)  ,  (At)2  (i-B)  v(t)  +  (At)20  V(t+At) 


(2.  19) 


(2.20) 


These  relations  are  used  to  update  the  motions  (u,  u,  u  *,  v,  v,  V)  for  each  time 
increment  At  at  ail  generic  points  in  the  space  domain  (t.e.,  mass  points  of  the 
model).  Assuming  the  motions  at  time  t  to  be  known,  the  process  of  updating 
consists  of  the  following  steps: 

(1)  At  time  t,  the  motions  of  all  mass  points  are  presumed  to  be 
known.  Compute  u(t+At)  and  v(t+At)  for  all  mass  points  from  the  appropriate 
equations  of  motions;  e.q.,  Eqs.  (2.1)  and  (2.2)  for  the  plane  strain  case,  using 
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the  stresses  at  time  t. 


(2)  Calculate  li(t^At),  v(t+At)  and  u(t+At),  v(t+At)  for  all  mass 
points  using  £qs.  (2.19)  and  (2.20),  respectively. 

(3)  Compute  new  strains  and  stresses  associated  with  the  velocities 
and  displacements  calculated  in  (2),  and  determine  u(t+&t)  and  v(t+At)  from  the 
appropriate  equations  of  motions  again. 

(4)  Compare  u  (t+At )  and  v(t+At)  of  steps  (3)  and  (1),  and  repeat  steps 
(2)  and  (3)  if  necessary;  otherwise,  increment  the  time  step  and  repeat  the 
process . 

2.3  Stability  Requirements  and  Analysis 

The  discretized  space  formulation  and  step-wise  numerical  integrator 
described  above  are  intended  to  yield  discrete  time  functions  of  the  relevant 
physical  quantities  at  all  generic  points  In  the  space  domain.  Unfortunately, 
not  every  numerical  solution  to  the  discretized  problem  is  necessarily  valid; 
there  are  certain  mathematical  requirements  that  must  be  satisfied  to  insure  a 
valid  solution.  In  other  words,  a  physical  discretization  of  the  space  and  time 
variables  is  not  sufficient  to  guarantee  a  valid  numerical  solution;  unless  certain 
mathematical  conditions  are  complied  with,  the  calculations  may  yield  completely 
erroneous  and  senseless  numerical  results.  This  is  in  contrast  to  many  statical 
problems  in  which  no  such  mathematical  conditions  are  required.  An  Important 
consideration  in  dynamical  calculations  is  the  assurance  of  stability  of  the 
numerical  scheme,  meaning  that  the  errors  associated  with  the  approximate 
numerical  solution  do  not  grow  with  time.  For  the  formulations  described  in 
Sects.  2.2  and  2.3,  the  required  stability  conditions  are  the  same  as  those  of 
finite  difference  schemes,  for  which  methods  of  stability  analyses  have  been 
well  developed  [  6,  8,  10,  12], 
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In  particular,  for  the  numerical  scheme  developed  herein  for  general 
problems  of  wave  propagation,  the  requisite  stability  conditions  can  be  obtained 
through  the  Fourier  transform  method  [12];  this  usually  yield  the  von  Neumann 
necessary  condition  for  stable  calculations.  A  complete  analysis  of  stability 
for  the  numerical  schemes  used  herein  can  be  found  In  Ref.  [II];  we  summarize 
herein  the  principal  results  for  plane  strain  and  axl-symmetrlc  media. 

As  it  turned  out,  the  stability  of  the  difference  schemes  depend  on  the 
value  of  0.  For  example,  in  the  plane  strain  case,  values  of  3  <  4  yields  only 
conditional  stability;  whereas,  for  9  >  i  the  resulting  difference  scheme  Is 
unconditionally  stable.  The  von  Neumann  condition  Is  only  a  necessary  condition 
for  stable  calculations;  however,  for  practical  purposes  this  is  often  sufficient. 

The  available  methods  of  stability  analysis,  including  the  Fourier 
transform  method,  are  based  on  an  explicit  set  of  finite  difference  equations. 

We  can  show  that  the  numerical  scheme  used  herein  Is  equivalent  to  an  explicit 
system  of  difference  equations.  Eq.  (2.20)  can  be  applied  to  time  t+2&t,  yielding 

u  (t  -2At)  =  u(t+At)  +  At-u(t+At)  +  (At)2(r-0)  tf(t+At)  +  (At  )20  ii(t+2At) 
an4  similarly  for  v(t+2At). 

Subtracting  this  from  Eq.  (2.20),  and  using  Eq.  (2.19)  to  eliminate  the  u  terms, 
we  obtain 

9  u' (t+2At )  +  ( 1 -2S )o  (t+At )  +  0  u(t) 

I 

-  - r  [  u(t+2At)  -  2u  (t  +At )  +  u(t)  ]  (2.21) 

(At)2 

Sim! larly, 

9  v  (t  +2  At )  +  (1-28)  v  (t+At )  +  8  v(t) 

I 

*  - r  [  v(t+2At)  -  2v  (t+At)  +  v(t)  ]  (2.22) 

(At)2 


18 


Thus,  the  Newmark  fj-integration  procedure  is  equivalent  to  solving  an  explicit 
system  of  difference  equations  with  the  time  discretization  indicated  in  Eq.  (2.21). 

Stability  of  a  numerical  solution  means  that  the  error  in  the  approxima¬ 
tion  remains  bounded  (does  not  grow  Indefinitely)  at  all  time  steps.  Since  the 
solution  vector  u(t)  at  any  time  t  of  the  resulting  system  of  difference  equations 
is  approximate,  it  will  contain  some  error,  say  error  vector  E.  It  can  easily  be 
shown  that  the  error  term  will  also  satisfy  the  homogeneous  part  of  Eqs.  (2.21)  and 
(2.22).  If  E  is  expanded  into  a  Fourier  series,  each  term  in  the  series  will 
individually  also  satisfy  the  homogeneous  system  of  difference  equations  if  the 
equations  are  linear.  A  typical  term  in  the  expansion  may  be  written  as 
exp  [  i  (ujjpAx  +  (ij^q'Ay)  ]  5°;  where  uuj  and  uij  are  frequencies,  p  and  q  are  each 
number  of  space  mesh,  n  is  the  number  of  time  steps,  and  $  is  the  modulus  of  the 
error.  For  the  difference  scheme  used  herein,  the  determinantal  equation  of  the 
homogeneous  system  yields  a  system  of  quadratic  equations  In  5-  The  error  term, 
therefore,  will  be  bounded  with  Increasing  time  steps  n  If  |{|  <  1.0;  hence, 
stability  of  the  solution  scheme  Is  assured  If 

lj|  <  1.0  (2.23) 


Stability  Conditions  for  Plane  Strain  Propagation  ---  For  elastic  wave 
propagation  under  plane  strain  conditions,  Eqs.  (2.21)  and  (2.22)  (after  using 
Eqs.  (2.6)  and  (2.7)  )  become 
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Th*  details  of  the  solution  of  Eqs.(2.26)  and  (2.27)  are  described  In  R»f.  (11]; 

the  resulting  solution  for  |^|  arc  shown  as  a  function  of  (c.*  In  fig.  2.1*. 

c  ay 

W,  observe  that  stability,  as  required  by  Eq.(2.29),  depends  on  the  value  of  0 
used  In  the  iicwrurk  integrator,  for  all  geometrically  meaningful  values  of  Bl 
l.r.,  f  «  0,  "jT~,  “ r— ,  and  -jj“,  the  corresponding  conditions  required  for 

stability  are  tabulated  below: 
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In  general,  the  stability  condition  is. 
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and,  unconditionally  stable  if  0  >  ^  . 

For  plane  strain  problems  formulated  with  the  non-Cartesian  model;  i.e., 
Eqs.(2.l3)  and  (2.14),  the  corresponding  stability  conditions  can  also  be  determined 
f  II]  through  the  above  Fourier  transform  method.  The  results  are: 
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whereas  for  S  £  ,  It  Is  unconditionally  stable. 


Stability  Conditions  for  Ax  I -symmetric  Propagation - For  wave 

propagation  under  axially  symmetric  condition,  Eqs.  (2.21)  and  (2.22), 
respectively,  yield 
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whe  re 


ana  lysis 
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wi  th, 


u.  ,  ’u(p  sr,  q  4z,  n  At) 

*  tJ 


wk,J  '  W*P  4r*  ^  4*.  n  At) 

On  the  basis  of  Eqs.  (2.28)  and  (2.29),  the  Fourier  method  of  stability 
leads  to  the  following  determlnantal  equations: 
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and  Cj  are  complex  quantities,  which  are 
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In  which  v  Is  the  Poisson  ratio. 
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For  generic  points  that  are  far  fron  the  axis  of  symmetry  (i.e.,  p  —  •  ), 

Cj  -•  Bj  and  —  B^,  and  Eqs.  (2.30)  and  (2.31)  approach  Eqs  (2.26)  and  (2. 27)  ,  respectively; 
hence,  the  corresponding  stability  conditions  tend  toward  those  for  the  plane  strain 
propagation. 

For  generic  points  that  are  close  to  the  axis  of  syrmetry  (I.e.,  finite 
p),  the  stability  analysis  requires  the  solution  of  Eqs.  (2.30)  and  (2.31)  for  |*  ] . 

The  solutions  are  given  by  ( 1 1  ]  . 

I  -  (i  -  »)C  ♦  */  0  -  (?  -  i)c12  -  (I  -  BC)2 

?«  -  (2  32) 

(I  *  0C) 

where  C  Is  Cj  or  Cj»  which  are  complex  quantities.  The  maximum  values  of  |»  | 
have  been  evaluated  numerically;  the  results  are  presented  in  Figs.  2.5  through  2.10 
which  show  the  envelopes  of  |{j  versus  K  for  all  geometrically  meaningful  values 
of  0  and  various  Aotsson  ratio  v.  An  exanl nat l on  of  these  figures  reveal  the 
fol lowing ; 

(1)  In  all  cases,  fg|  >  1.0;  this  means  that  the  general  numerical 
scheme  Is  locally  unstable  when  applied  to  axl-symmctrlc  wave  propagation. 

(2)  As  i  Increases,  the  ’degree"  of  Instability  decreases  (I.e.,  |f| 

Is  not  much  greater  than  1.0),  but  In  no  case  can  stability  be  assured  In 
accordance  with  Eq.  (2  .23). 

(3)  As  v  Increases,  the  degree  of  Instability  worsens. 

(U)  As  p  Increases,  the  modulus  |*|  tends  toward  the  value  1.0. 
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Stability  of  the  original  scheme,  therefore,  cannot  be  assured,  on 
the  basis  of  Eq .  (23),  for  ax i -symmet r ic  calculations.  However,  conditional 
stability  of  the  difference  scheme  can  be  achieved  through  the  Introduction  of 
a  dissipative  mechanism,  in  the  form  of  artificial  viscosity  [12),  into  the  original 
system  of  equations,  such  that  the  equivalent  system  of  differential  equations 
become : 


*  (crr  +  qr)  d 


ar 


rz 
dz 


P  u 


(2.33) 


rz 


*r 


» (ffz  +  qz) 


dz 


rz 

+  —  =<  p  a 
r 


(2.3*0 


where  qr  and  q^  are  pseudo-viscous  stresses  In  the  r  and  z-dl rectlons,  respectively, 
corresponding  to  an  artificial  viscosity  coefficient  F,  given  as  follows: 
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Clearly,  without  the  terms  7—  and  t —  ,  the  discrete  forms  of  Eqs.  (2.33)  and 

0“  oz 

(2.34)  would  reduce,  respectively,  to  Eqs.  (2.15)  and  (2.16).  With  the  Introduction 

of  artificial  viscosity,  discretized  forms  of  ^r  and  must  be  added  to 

•r  8z 

Eqs.  (2.15)  and  (2.16) .  This  may  be  done  by  applying  backward  difference  to  iju  and 

ir 

of  Eqs. (2.35)  and  (2.36);  thus, 
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*  The  corresponding  forward  difference  will  lead  to  an  unstable  scheme;  hence, 
should  be  avoided  [It]. 
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and, 
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Then,  applying  centered  difference  for  the  space  derivatives,  we  obtain  for  the 
artificial  viscosity  terms, 
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Eqs.  (2 .35a) and  (2.36a)  should  be  added,  respectively,  to  Eqs.  (2.28)  and  (2.29).  The 
Fourier  transform  method  applied  to  the  resulting  systems  of  equations  then  yield 
the  following  determlnanta 1  equations  for  the  modified  difference  scheme: 

».?3  tbi?2  +cl?  *«l  -°  (2-37> 

a2?3  +  b252  +  c2S  +  d2  =  0  (2.  38) 

whe re : 

a  |  =*  1  +  3Cj  +  j  a2  +  1  +  0C2  +  BAj 

b,  -  -2  +  (I  -  20 ) C |  +  (I  -  3B)A,  ;  b2  -  -2  +  (1  -  23)C2  +  (1  -  3B)a2 

Cj  -  1  +  0Cj  -  (I  -  30 ) A |  ;  c2  =  1  +  0C2  -  (1  -  3S)A2 

dj  =  — 0A |  ’  d2  =  "®A2 

in  which  Cj  and  C2  are  defined  following  Eqs.  (2.30)  and  (2.31),  and 
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The  roots  f  of  Eqs.  (2.37)  and  (2.38)  can  only  be  determined  numerically 
the  results  for  various  values  of  T  and  v,  and  at  varying  distances  from  the  axis 
of  symmetry  (l.e.,  different  p),  are  summarized  in  Ftgs.  2.11  through  2.16 
Only  the  maximum  moduli  of  5  are  given  in  these  figures.  Examination  of  the  plots 
shown  in  these  figures  indicate  that  depending  on  the  values  of  0  and  p,  there 
is  a  range  of  values  of'  F  for  which  the  von  Neumann  condition  j^l  <  1.0  is 
satisfied  for  values  of  v  <  0.50:  and  hence  stability  can  be  conditionally 
assured.  However,  for  v  =  0.50  (Incompressible  medium)  the  use  of  artificial 
viscosity  does  not  lead  to  a  stable  scheme;  see  Figs.  2.13  through  2.16.  Figs. 

2.11  through  2.15  also  show  that  as  p  or  8  Increases,  the  range  of  stability 
corresponding  to  a  given  value  of  r  Improves. 

From  these  figures,  It  can  be  seen  that  for  specified  values  of  p,  8, 
and  v t  the  stable  range  of  K  (l.e.,  range  of  K  for  which  |5|  <  1.0)  varies  with 
T.  Also,  corresponding  to  each  value  of  T,  there  Is  a  nontrivial  value  of  K 
(l.e.  t  0)  for  which  |r|  a  1.0;  this  value  of  K  can  be  called  Kcr»  Values  of 
Kc(_  are  plotted  versus  T  in  Figs.  2.17  through  2.22;  It  can  be  observed  from  these 
figures  that  depending  on  p,  0,  and  v,  K  attains  a  maximum  value  for  certain  T. 

c  “  i 

This  value  of  F  represents  an  optimal  value  of  artificial  viscosity,  and  is 
recommended  for  practical  applications. 
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III.  DEVELOPMENT  OF  TRANSMITTING  BOUNDARY 


3 . I  Need  for  Transmitting  Boundary 

The  discrete  method  of  formulation  described  in  Chapter  II  is  a  general 
numerical  method  for  determining  the  physical  quantities  at  the  discrete  generic 
points  in  space  as  well  as  In  discrete  instants  of  time;  in  other  words,  a 
solution  obtained  through  such  a  formulation  is  a  discrete-variable  function  in 
space  and  time.  Such  a  discrete-variable  solution  is  necessarily  limited  to  a 
finite  region  of  the  physical  space;  a  solution  for  a  space-time  problem  then 
consists  of  the  relevant  quantities  within  a  finite  space-time  domain. 

The  space  domain,  therefore,  must  always  be  bounded  or  terminated  by 
suitable  boundary  conditions.  If  these  boundaries  are  well  defined,  such  as  a 
regular  stress  or  displacement  boundary,  then  the  problem  Is  completely  described 
in  the  discrete  sense.  However,  If  a  problem  involves  an  infinite  or  semi-infinite 
medium,  then  it  is  not  possible  to  describe  the  solution  over  the  entire  space; 
to  do  so  would  require  an  infinite  number  of  discrete  points.  The  alternative  Is 
to  terminate  the  medium  at  some  appropriate  location  with  an  artificial  boundary 
that  will  reproduce  essentially  the  same  effect  as  that  of  the  infinite  medium 
beyond  It.  In  other  words,  the  artificial  boundary  must  be  such  that  all  incident 
waves  are  transmitted  through  it  without  reflection,  as  If  the  material  were 
continuous  and  the  boundary  were  not  present.  Such  a  termination  may  be  called 
a  "transmitting  boundary"  or  a  "nonreflecting  boundary". 

Alternatively,  a  ca leu latlona 1  scheme  may  be  developed  to  serve  the 
same  purpose  as  that  of  an  explicit  transmitting  boundary  by  taking  Into 
consideration  the  theoretical  effects  of  an  Infinite  or  semi-infinite  region. 

In  practice,  transmitting  boundaries  are  often  required  for  reasons  of 
economy  in  calculations.  For  Instance,  in  the  practical  prediction  of  ground 
motions  induced  by  nuclear  bursts,  such  transmitting  boundaries  are  needed  in  a 
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number  of  situations,  including  the  following: 

(1)  Calculation  of  ground  motions  in  the  distant  regions  of  a  burst 
where  generally  outrunning  conditions  prevail. 

(2)  Calculation  of  complex  structure -medium  interaction  under  a  ground- 
transmitted  disturbance  and/or  air-blast  forces  resulting  from  a  nuclear  burst. 

in  the  first  instance,  the  motions  and  stresses  at  relatively  shallow  depths  in 
tie  distant  outrunning  regions  arc  of  interest;  unless  the  calculations  of  the 
motions  originating  from  the  source  are  limitad  to  the  upper  surface  strata, 
terminated  properly  with  a  suitable  nonreflecting  boundary,  the  calculations  must 
be  extended  to  the  regions  at  great  depths  In  order  to  avoid  the  artificial 
reflection  at  the  bottom  from  reaching  the  outrunning  regions.  In  the  second  case, 
a  transmitting  boundary  will  permit  the  isolation  of  the  region  around  the  structure 
and  the  calculations  can  be  confined  within  this  region  without  the  unwanted 
reflections  from  the  terminating  boundaries  of  tho  region;  unless  this  can  be  done, 
the  calculations  of  structure -medium  interaction  effects  must  necessarily  be  ■ 
performed  with  very  coarse  grid. 

3 . 2  Previous  Work 

*.v 

One  of  the  first'  solutions  Involving  the  use  of  a  transmitting  boundary 
was  performed  for  a  problem  involving  outrunning  ground  motions  [3 ] .  Tha 
proh'.  •  vo  s  the  detennination  of  ground^notion  histories  due  to  a  surface 

nuclear  bv.  a  layered  system.  The  problem  was  modeled  to  consist  of  three 

layers  of  cl*  ~>c  material,  with  the  last  layer  assumed  to  have  infinite  depth. 

The  input  consists  of  the  direct  ground  shock  effect  applied  in  a  crater  plus 
the  expanding  air-blast  pressures  applied  at  the  surface. 


*  To  the  knowledge  of  the  authors,  this  was  the  first  successful  use  of  a 
transmitting  boundary  in  numerical  calculation  of  wave  motions. 
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In  order  to  obtain  reasonably  reliable  results  In  the  outrunning  regions 
(at  large  distances  from  ground  aero),  the  depth  of  the  spatial  region  Involved  In 
the  calculations  must  be  limited,  otherwise  the  amount  of  storage  and  calculation 
*.lme  would  be  excessive;  to  do  this,  the  bottom  was  terminated  with  a  boundary  that 
simulated  the  Infinite  depth.  An  early  version  of  a  transmitting  boundary  was 
employed  for  these  calculations  [3].  This  Is  the  forerunner  of  the  trans* 
mlttlng  boundary  subsequently  refined  and  f  rther  developed  In  the  studies 
described  herein. 


3.3  Theoretical  Basis  of  Trensml tt tng  Boundary 

The  concept  of  a  transmitting  boundary  or  a  'Vtonref lectlng  boundary"  for 
wave  propagation  calculation  as  conceived  In  this  project  can  be  described  first 
(for  the  sake  of  clarity)  for  the  elastic  plane  one-d (mens Iona  1  case.  This  same 
concept  can  be  extended  to  higher  space  dimensions  and  will  be  described  subsequently. 

Consider  a  seml-Inf Inlte  bar  subjected  to  a  plane  stress  pulse  P(t)  from 
the  free  end,  as  shown  In  Fig.  3.1a.  The  discrete  Idealisation  of  the  bar,  based 
on  the  same  space  discretization  as  those  for  the  plane  strain  and  axlsymmetrlc 
media.  Is  shown  In  Fig.  3.1b.  We  observe  that  the  equation  of  motion  of  a  general 
Interior  mass  point  I  at  time  t  Is, 


*‘(1  +  I)  -  0*0  -  I)  t 

-  -  0ff*(J) 

A* 


(3.D 


The  corresponding  equation  for  the  mess  point  b  at  the  terminating  boundary  Is 
slml larly, 

o*  (b  +  I )  -  o*  (b  -  I )  „ 

-  -  p  u*(b)  (3.2) 

Ax 

However,  for  a  mass  point  on  the  boundary,  the  stress  a*  (b  ♦  I)  Is  not  known  and 
cannot  be  computed  In  the  usual  manner.  But  If  the  velocity  of  propagation  Is  c, 
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It  c.in  be  observed  from  Fig.  3.1b  that  the  travel  tlM  of  a  strata  wave  fra  ona 
stress  point  to  the  next  is  exactly  equal  to  ona  increment  of  transit  tlaa  h  * 

A*/ c.  Hence,  it  the  stress  at  stress  point  (b-l)  at  tlaa  (t-h)  Is  ot"**(b-l), 
then  this  stress  must  be  the  stress  at  point  (b+t)  at  tlaa  t;  thus* 

-*(b  •  I)  -  rXmh(b  -  I)  0.1) 

Therefore,  In  order  to  simulate  the  Infinite  medium  the  aquation  of  ant  I  on  of  the 
boundary  mass  point,  Eq.  (3.2),  must  be, 

ffl"^b*  I)  -  (b  -  I)  t 

- -  p  st(b)  (l.%) 

As 

We  might  emphasize  that  for  plana  elastic  propagation,  Eq.  (J.4)  will  slaw  I  ate 
an  Infinite  medium  exact  Iv  (In  the  sense  that  there  Is  no  additional  approx  law  t  ton 
other  than  those  of  the  centered  finite  difference).  However,  the  extension  to 
higher  space  Is  not  obvious  on  this  basis.  For  this  latter  purpose',  we  shell  derive 
Eq.  (}.4)  from  another  consideration. 

In  discrete  steps,  the  transmission  of  the  stress  wave  may  be  considered 
as  the  transfer  of  the  D'AleaAert  force  at  t  mess  point  from  one  stress  point  to 
the  next  as  time  Increases  In  Increments  of  h.  From  this  standpoint,  the  D'AleaAert 
force  on  the  mass  point  I  at  t law  (t-h)  Is  equal  to 

-  A  (  f+il  *  I}  -  d*’h(l  -  I)  J 

where  A  Is  the  area  of  the  bar.  Transmission  of  this  force  will  give  rise  to  a 
change  In  the  force  of  stress  point  (I  ♦  I)  from  time  (t  -  h)  to  tlaw  t  equal  to 
the  D'Alembert  force  on  mass  point  I  at  time  (t-h);  l.e.. 
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Therefore 


(3.3) 


A  [  <,*(}  +  I)  -  <jt_h0  +  I)  ]  »  -A  (  a*-h(l  +  1)  -  <Tt_h(i  -  I)  ] 

.  <7*  (I  +  l)  -  at_h(i  -  D 


Thus,  for  the  boundary  mass  point,  i  *  b,  Eq.  (3.4)  is  again  obtained  from  Eq. 

(3.1)  using  Eq.  (3.3). 

3 .4  Transmitting  Boundary  in  Plane  Strain 

The  extension  of  the  above  concept  of  the  transmission  of  D'Alembert 
forces  to  higher  space  will  now  be  illustrated  for  the  two-dimensional  space  under 
plane  strain  condition.  Consider  the  lumped-parameter  idealization  of  a  two-space 
as  shown  in  Fig.  3.2.  Assume  that  the  region  is  terminated  at  j  *  b.  The  equations 
of  motion  of  the  mass  point  (i,b)  are: 

oJ(l  +  1  >b)  -  o*(l  -  1,b)  T*y(»,b  +  1)  -  fjy(i,b  -  1) 

Lx  Ay 

-  p  ii*  (I  ,b)  (3.5) 


and. 


oJ(l,b  +  I)  -  oJ(t,b  -  I) 

. -  ■  ■  ■  — — ■  + 

Ay 


Ax 


*  p  v*  (I  ,b) 


(3.6) 


From  Fig.  3.2,  we  observe  that  the  following  physical  quantities  are  not  defined 
in  the  above  equations  of  motions,  and  cannot  be  determined  in  the  usual  manner: 
(I )  The  stresses  g  (I ,b  +  I )  and  t  (I ,  b  +  1 )  for  all  i  ;  and 

y  xy 

(ii)  The  strain  •  (i  +  l,b)  and  t  (i  +  l,b)  for  all  1.  <  (i  +  l,b) 

y  xy  y 

is  required  to  determine  0^(1  +  l,b). 


33 


These  quantities  can  be  determined  from  a  consideration  of  the  transmission  of 
the  D'Alembert  forces  on  the  boundary  mass  points. 

Consider  first  the  mass  point  (l,b).  If  the  stresses  are  assumed  to 
propagate  at  the  dllatatlonal  velocity  of  the  medium  c^  ,  then  from  the  trans¬ 
mission  of  the  D'Alembert  force  on  this  mass  point  in  the  x-directlon,  we  have 


-  {  Ay  I  <r*“h(l  *  MO  -  o*“h(l-1,b)  ] 


-  Ax  [  T*-h(l,b  ♦  l)  -  T*'h(l,b  -  I)  ]  } 


*  Ax  [  t*  (l.b  +  I) 

*7 


t-h 

Txy  (t,b  +  I)  ] 


from  which  we  obtain, 

T*y,(l,b  ♦  I)  f  cj‘h(l  ♦  l,b>  -  oJ"h(l  -  l,b)  ] 

Ax 

♦  TiCh(,'b  *  0  (3.7) 

Ay 

where  h  ■  —  ;  the  unit  dllatatlonal  transit  time. 

ed 

Similarly,  consideration  of  the  D'Alembert  force  In  the  y-direction  yields, 

-  {  Ax  [  *y"h(l#b  ♦  |)  -  <Ty"h(l ,b  -  I)  ] 

-  Ay  [  rJjNl  ♦  l,b)  -  T^h(t  -  l,b>  ]  } 

•  Ax  {  *J(l,b  ♦  I)  -  oJ’h(l,b  +  I)  ] 

from  which, 

oi,(l,b  ♦  I)  -  m*’h(l,b  -!)♦—[  T*’h(l  -  l,b)  -  r*“h(l  +  I»b)  )  (3.8) 

7  1  7  Aw  *7  *7 
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Based  on  the  same  assumption  that  stresses  are  propagating  at  the 


dllatationa!  velocity  c^,  the  D'Alembert  force  In  the  y-dlrectlon  of  mass  point 
(I  +  I  ,b  -  I)  at  time  (t-h)  must  be  equal  to  the  change  In  the  vertical  force  at 
stress  point  (I  +  l,b)  from  time  (t-h)  to  t;  thus, 

-  {  Ax  [  crj"h(l  +  l,b)  -  0J’h(l  +  l,b  -  2)  ] 

+  Ay  I  T*‘h(l  +  2,b  -  I)  -T*‘h(l,b  -  l)  ]  } 

-  AX  [  oj  (I  +  l,b)  -  oJ‘h(l  +  l,b)  ] 

Therefore, 

ex*  (I  +  l,b)  -  o*‘h(!  +  l,b  -  2)  -  —  [  T*‘h(l  +  2,b  -  I)  -  r‘‘h(l,b  -  I)]  (3.9) 

y  ■  y  ^  *y  *y 

Knowing  <y*(l  +  l,b),  the  strain  «y(l  +  l,b)  Is  obtained  from  the  stress-strain 
equation.  This  will  then  permit  the  calculation  of  the  stress  c*j(l  +  l,b). 

<7*1  (I  -  I ,b)  Is  found  similarly  through  the  transmission  of  the  vertical  D'Alembert 
force  of  mass  point  (I  -I,  b  -  I). 

Considering  the  D'Alembert  force  of  mass  point  (I  +  l,b  -  I)  In  the 
x-dlrectton,  we  obtain 

-  {  Ay  l  o‘"h(l  +  2,b  -  I)  -  -  I)  J 

+  Ax  [  T*’h(l  +  t,b)  -  T*'h(i  *  l,b  -  2)  ] 

*y  *y 

-  AX  I  T*y (I  ♦  l,b)  -  T*“h<l  ♦  l,b)  J 


from  which. 
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(3.10) 


%vlO‘».b)  =  *l,b-2) - I  ~'"h(l*2,b-l)  -  ®*‘h(l,b-l)  ) 

xyi  xy  Ax  * 

Sint larly , 

t'  ,(M,b)  -  r*'h(i*l,b-2) - (  ***h(l,b-l)  -  o‘*h(l-2,b-l)  ]  (3.11a) 

y  V  Ax 

The  stresses  given  In  Eqs.  (3.7)  through  (3.1V)*  which  are  required  In 
Eqs.  (3.5)  and  (3.6),  would  be  the  correct  stresses  If  the  velocity  of  propagation 
of  these  stresses  were  Indeed  cd<  However,  In  a  two-dimensional  solid  mad  tun,  It 
Is  known  that  certain  stresses  propagate  at  the  shear  velocity  c%.  Alternatively* 

If  the  stresses  through  the  boundary  are  assumed  to  propagate  at  the  shear 
velocity  cs  ,  consideration  of  the  corresponding  D'Aleafcert  forces  at  the  boundary 
mass  points  would  yield  the  following  stresses*  Instead  of  Eqs.  (3.7)  through  (3. II): 


t 

Txy2 

(1  ,b  +1 )  - 

*Tt-k 

xy 

(»,b- 

AV  ,-k 

1) - [  <J*  k(l+l,b)  - 

Ax 

o*"k(l-l,b)  ) 

(3.7e) 

*J20.b*l)  -  <rj’k(l,b 

-1)  ♦ 

—  i  %ck(M»b)  -  J 

Ax  y  7 

(3.8a) 

+l,b)  »«-y"k(l+l 

,b-2) 

Ay  r-k 

-  —  (  T*;k(l*2,b-I)  - 
Ax  ** 

T*"k(l,b-I)  ) 

(3.90 

t 

Txy2 

(1  *l,b)  - 

rt-k(l* 

xy 

1  ,b-2 

Ay  f.k 

) - (  **  k (1  *2,  b-l) 

Ax 

-  o*"k(l*b-l)  ) 

(3.10a) 

t 

Txy2 

(l-l,b)  - 

T*’k(l- 

xy 

l,b-2) - (  p*"k(l,b-l)  -  ffj’k(l-2,b-l)  ] 

Ax 

(3.11a) 

which 

k  -  Jfc 

r 

,  the 

unit 

shear  transit  time. 

However,  the  transmission  of  the  D'Alembert  forces  will  not  be 

exclusively  at  the  velocity  c  nor  exclusively  at  c  ;  Invariably*  It  will  be 

d  * 
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a  combination  of  the  two  propagation  velocities.  Accordingly,  the  correct 
stresses  to  be  used  in  Eqs.  (3.5)  and  (3.6)  should  be  some  combinations  of  Eqs. 
(3.7)  through  (3.11)  with  Eqs.  (3.7a)  through  (3.11a).  Unfortunately,  there  is 
no  theoretical  basis  for  combining  these  stresses.  Through  various  computational 
experiments,  it  was  determined  that  the  following  combination  is  quite  suitable 
and  sufficiently  accurate  for  practical  purposes. 


For  a  direct  stress. 


and  for  a  shear  stress. 


3  •  +  °  -  •>  'U 


Txy  -<'-•)  Txy I  +  •  Txy2 


(3.12) 


(3.13) 


In  which  a  Is  the  fraction  of  the  total  O'Alembert  force  contributed  by  the  direct 
stresses  when  the  velocity  of  propagation  is  assumed  to  be  c.  ,  whereas  B  is  the 

d 

fraction  of  the  shear  stress  contribution  to  the  total  D'Alembert  force  when  the 
propagation  velocity  is  assumed  to  be  cs  .  For  example,  when  applying  Eqs.  (3.12) 
and  (3.13)  to  the  stresses  at  (t,b+l), 


I  aj‘h(t,b+l)  -  «rj"h(l,b-l)  { 

|  <7y"h(i,b+l)  -  <rj‘h(l,b-l)  |  |  r^h(i+l,b)  -  TJ‘h(l-l,b)  | 

and, 

I  rj"k(l,b+l)  -  TJykO,b-l)  | 

!  rjyk(i,b+l)  -  T*’k(i,b-I)  |  +  |  o*-k(l*l,b)  -  <y*'k(i-l,b)  | 


(3.14) 


(3.15) 


We  might  observe  that  for  a  plane  dllatational  propagation,  Eqs.  (3.14) 
and  (3.15)  become 


or  »  1 .0  and  9-0 
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and  tf  the  propagation  is  also  in  a  direction  normal  to  the  transmitting  boundary, 
the  transmitted  stresses  become 


and, 


T*y(i,b+I)  *  0 

ey(l,b+l)  -  oyl(?,b+1) 

-  CTy“h(l,b-l) 


which  is  the  same  as  Eq .  (3.3)  for  one-dimensional  plane  propagation. 

On  the  other  hand,  for  plane  shear  propagation,  Eqs,  (3. 1*0  and  (3.15) 

became 


or  »  0  and  9  ■  1 .0 


and  if  the  direction  of  propagation  is  normal  to  the  boundary,  the  transmitted 
shear  stresses  are 

%y(,’b*n  ■Txy2-^k(,'b-') 
and  cry  (I  ,b+l )  ■  0 

which  are  the  results  for  the  one-dlmens Iona  I  plane  shear  propagation. 

Special  Boundary  for  Linearly  Elastic  Media  —  The  transmitting  boundary 
described  above  is  applicable  for  problems  In  which  inelastic  behavior  may  occur 
anywhere  within  the  region  terminated  by  these  transmitting  boundaries.  However, 
for  problems  involving  purely  elastic  behavior,  the  lumped-parameter  model  of  Fig. 
3.2  can  be  simplified  to  that  of  Fig.  3.3,  in  which  each  mass  point  defines  either 
the  vertical  or  horizontal  motions  of  the  model  only;  accordingly,  the  direct  and 
shear  strains  and  stresses  are  also  defined  only  at  alternate  stress  points,  as 
shown  in  Fig.  3.3. 
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Consistent  with  the  lumped-parameter  model  of  Fig.  3.3,  the 
transmitting  boundary  can  also  be  simplified.  Specifically,  In  this  case,  we 
see  from  Fig.  3.3  that  the  vertical  D'Alembert  force  of  mass  point  (l,b) 
yields  Eqs.  (3*8)  and  (3.8a);  whereas,  the  horizontal  D'Alembert  force  of  mass 
point  (1+1, b-1)  yields  Eqs.  (3.10)  and  (3.10a).  On  the  basis  of  Eqs.  (3.12) 
through  (3 . 1 5 ) >  the  stresses  required  In  the  equations  of  motions  of  the 
boundary  mass  points  are,  therefore,  completely  defined. 

3 .5  Transmitting  Boundary  In  Ax  I -symmetric  Condition 

The  concept  of  the  step-wise  transmission  of  D'Alembert  forces  can  be 
readily  extended  to  the  development  of  an  axi-symmetric  transmitting  boundary. 

Consider  the  lumped-parameter  model  of  an  axl-symmetrlc  half-space 
shown  In  Fig.  3.4,  which  is  terminated  at  J  »  b.  We  observe  that  the  motions 
of  all  mass  points  at  J  *  b  and  j  -  b-1  cannot  be  determined  in  the  usual  manner; 
that  Is,  the  equations  of  motions  of  these  mass  points  will  contain  stress  (or 
strain)  components  which  cannot  be  defined  In  the  usual  way  In  terms  of  differences 
In  displacements.  For  example,  the  equations  of  motions  of  mass  point  (l,b)  at 
time  t  are: 


<Jr(1+1,b)  -  err (1-1  ,b)  Trz(l,b+!)  -  Trz(l,b-1) 


tr 


Az 

ffr(!  ,b)  -  (I  ,b) 

r(») 


P  u  (I  ,b) 


and  , 


rrz0+1,b)  -  rr2(l-1,b)  trz(1,b+l)  -  cz(l,b-l)  Tpz  (I  ,b) 


Ar 


Az 


r(D 


*  o  w(l,b) 


(3.16) 


(3.17) 


The  following  stresses  cannot  be  determined  In  the  usual  manner: 
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In  Eq  .  (3.16)'. 

r  (l,b’l)  is  completely  unspecified;  whereas,  i?r(i+l,b)  and  »r(l-l,b) 
cannot  be  cimDuted  because  e^(i*l,b)  and  «^(l-l,b)  are  not  defined. 

In  Eq.  (3.17); 

<7  (i.b-l)  is  completely  unspecified;  whereas,  r^fi'-ljb)  and  t  (1  —  1  ,b ) 
cannot  be  computed  because  is  not  defined  for  the  shear  strains  at  these 
stress  points. 

However,  these  stresses  can  be  determined  on  the  basis  of  the 
transmission  of  D'Alembert  forces  on  the  boundary  mass  points;  i  e. ,  the  mass 
points  at  j  »  b  and  j  *  b-l.  Specifically,  consideration  of  the  D'Alembert  force 
in  the  z-dlrection  of  mass  point  (1*1, b-l),  assuming  the  transmission  at  the 
veloc I ty  cd,  yields 

cr*.  (I  *-1  ,b)  -  o*‘h(l+l,b-2)  -  —  [  T*"h(l *2,b-l )  -  tj"h(l,b-1)  ] 
zl  z  Ar 

- —  T*'hO+l,b-l)  (3-18) 

r (I  *1 )  fZ 

Az 

where  h  »  —  . 

cd 

Knowing  this  stress,  the  corresponding  strain  «*(l+l,b)  can  be  determined  from 
the  pertinent  stress-strain  equation.  This  will  then  permit  the  calculation  of 
e*|(l‘l,b).  On  the  same  basis,  ff*|(l-l»b)  Is  determined. 

Consideration  of  the  D'Alembert  force  In  the  r-dlrectlon  yields. 


Az 


t*  .  (I  *1  ,b)  -  rJ‘h(l+l,b-2)  -  —  (  ffj'h(«*2,b-1)  -  «*”n(1 , b-l )  ) 


t-h  i 


Ar 


♦  A-2—  (  ff!*h(l*i,b-l)  -  c*‘h(l*1,b-l)  ] 

r  (1*1)  9  r 


(3.19) 
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Applying  the  transmission  (assuming  a  velocity  c^)  of  D'Alembert  forces 
to  mass  point  (i,b)  then  yields, 


^,(l,b+l)  =  o‘‘h(i,b-1)  -  —  [  T^h(i+I,b)  -  Tj"h(l-l,b)  ] 

Ar 


Az 


r  (i ) 


T^h(i,b) 


(3.20) 


and, 


TfzI(l,b+l) 


Tj:h(l,b-1)  -  —  [  <?J"h(l  -l,b)  -  a*“h(l-l,b) 

Ar 


+  —  [  ^’h0,b)  -  o*‘h(l,b)  ] 

r  (I )  9  r 


(3.21) 


Eqs.  (3.18)  through  (3.21)  were  obtained  assuming  that  the  propagation 
velocity  Is  c^  .  Alternatively,  If  the  velocity  of  propagation  Is  assumed  to  be 
cg  ,  the  corresponding  stresses  would  be 

<rz2(l+l,b)  ,  Tjz2(l+I,b)  ,  0*2  (I  ,b+l  )#  and  tJz2  (I  ,b+l )  , 

which  are  obtained  by  replacing  h  by  k  ■  Az/cs  in  Eqs.  (3.18)  through  (3.21). 

From  cr£2(!+l,b),  the  corresponding  strain  e*(l+l,b)  is  obtained,  from  which 
oz2(l+l,b)  Is  determined.  o^2(l-l,b)  is  obtained  on  a  similar  basts. 

The  correct  stresses  to  use  in  Eqs.  (3.16)  and  (3.17)  will  be  some 
combination  of  the  stresses  computed  on  the  basis  of  velocity  crf  with  the  stresses 
corresponding  to  velocity  cs  . 

Again,  on  the  basis  of  various  computational  experiments,  the  combinations 
suggested  in  Eqs.  (3.12)  and  (3.13)  were  found  to  be  also  suitable  for  axi- 
symmetric  calculations;  specifically,  in  this  case  we  have. 
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zl 


(I  -  y)  a 


z  2 


(3-22) 


and . 


t 

rz 


=  (I  -  B)  t* 


rz  I 


S  T 


rz2 


(3.23) 


whtT..  the  factors  y  and  ?,  when  applied  to  the  stresses  at  (i,b>l),  for  example, 
are  as  follows: 


|'/'h(i,b*l)  -  -r‘'h(i,b-l)| 


i"7  (*»b  +  l)  -  a*  (i,b-l)|  +  |  t^*zh(l '■I  ,b)-t^zh(i-1  ,b).  -  7Jf)TpZh  (' >b )  I 

...  (3.24) 


and , 


lT«k(,»b*l)  *  V"k0,b-l)| 


A  = 


l\’k  ('  »b *1  )-rJ-k  (i  ,b-l )  |  *  |— [9J_k  (I  H  ,b)-^_k  (i-l,b)]  +  777y(»J'k  (I  ,b)-<yj'k  (i 

Ar 

...  (3.25) 

We  can  observe  also  In  this  case  that  if  the  propagation  is  purely 
d  i  la  tat  i  ona  I  ,  Eqs.  (3.24)  and  (3-25)  become 


y  =  1.0  and  9=0 

and  if  the  propagation  is  normal  to  the  transmitting  boundary,  Eqs.  (3.18)  through 
(3.21)  would  yield  the 


and. 


Trz  (i  »b  *■• )  =  0 

<J*(?,bM)  =  n*|(i,b+l) 


=  o‘”h(|,b-l) 

i. 


>b)]| 
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which  are  the  results  expected  in  a  one-dimensiona I  plane  propagation. 

Special  Boundary  for  Ax ?- symmetric  Elastic  Media  ---  In  general,  two 
sections  of  mass  points  and  stress  points  arc  required  to  describe  an  ax i -symmetric 
solid;  the  grids  between  the  two  sections  are  shifted  by  half  a  nesn  length  in 
each  direction,  as  indicated  previously  and  shown  in  Fiq.  2.3. 

However,  for  linearly  elastic  media,  a  single  grid,  on  which  the  radial 
and  vertical  motions  are  defined  at  alternate  mass  points,  will  suffice  If  certain 
quantities  are  approximated  with  the  corresponding  averages  of  the  neighboring  mass 
points  (or  stress  points  as  appropriate).  Such  a  simplification  is  similar  to  that 
described  in  Sect.  3  .b  for  linearly  clastic  material  In  plane  strain. 

3-6  Transmitting  8oundarv  for  Elastic-Plastic  Propagation 

for  elastic-plastic  material,  the  transmitting  boundary  described  above 
can  be  used  without  modification  If  the  necessary  calculations  are  terminated  In 
a  region  where  the  material  can  be  assumed  to  be  clastic.  In  many  cases,  this 
should  be  sufficient  and  whenever  applicable  should  be  used. 

However,  the  above  transmitting  boundary  can  also  be  extended  to  Include 
clastic-plastic  propagation.  Basically,  this  extension  for  elastic-plastic  media 
must  Include  the  fact  that  there  are  more  than  one  dllatatlonal  propagation 
velocity;  In  particular,  the  elastic  component  of  a  wave  will  propagate  ahead  of 
the  plastic  component.  Also,  the  unloading  wave  may  propagate  at  Its  own  velocity 
(elastic  unloading  Is  assumed  In  this  study). 

Several  schemes  for  Including  the  effects  of  elastic-plastic  behavior 
were  investigated  In  this  study.  On  the  basis  of  these  Investigations,  the 
following  extension  or  modification  to  the  algorithms  described  In  Sec^.  3 •  ^ 
and  3*5  I*  suggested: 


Consider  the  plane-strain  case  (the  extension  to  the  axi-symmetrlc 
case  can  be  similarly  accomplished).  Referring  to  Fig.  3.2  ,  we  determine 
the  stresses  at  point  (i,b*l)  as  follows: 


Let  Y  be  an  indicator  at 


time  (t-hp),  such  that 


Y 


f  0  ,  if  the  material  Is  elastic 

i  1  ,  If  the  material  is  plastic 

-I  ,  If  the  material  Is  unloading 


Then  the  stress  s*(l,bl)  Is  determined  as  follows: 


.-J(l,b+I)  »  (l-X,)  (  ,v  -J'h(i,b-I)  +  (I-*)  crj'k(l,b-l)  ] 

+  X,  (  *  nt"h(l,b-l)  ♦  (!-»)  o*"k(l,b-») 

e 

*  (  ff*”hp(l»b-l)  -  ^*hp(l,b-l)  ] 

♦  —  -f  (l-X,)  [  *  r*‘h(l-l,b)  ♦  (l-»)  r*’k(l-l,b)  ] 

&X  1  2  xy  *y 

+  x2  [  ®  Tx^(,’'*b)  *  0_ty)  Txye(*“,,b) 
b)  - 

.  —  .  (I-K.)  [  #  Nl*l.b)  *  (1~.)  .Xy  (t-1  .b)  1 

A* 

+  X3  [  w  Txy^(,-'»b)  *  ('*»)  \£(,’,»b) 


(3.26) 
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Jn  which  h  =*  —  ,  the  ful  ly  plastic  unit  transit  time  with  c  -  /  --/r.  > —  : 

p  cp  P  V3(I-2vTo 

|  y  |  *♦  y 

and  X  -  j -  ;  where  the  subscripts  1,  2,  3  on  X  refer  to  the  stress  points 

(l»b-l),  (i-l,b),  and  (1+1, b),  respectively,  at  time  (t-hp). 

We  should  emphasize  that  If  the  material  remains  purely  elastic  (l.e., 
V  •  0),  Eq.  (3.26)  reduce  to  Eq .  (3.12);  whereas.  If  a  stress  point  is  elastic- 
plastic  (l.e.,  Y  ■  I),  the  elastic  precursor  will  be  transmitted  through  the 
boundary  at  the  elastic  velocity  and  the  plastic  component  will  be  transmitted 
at  the  fully  plastic  velocity  cp. 


Slml larly. 


d,b+i)  -  (i-x, )  [e 


+  x|  I  9  w(l,b-l)  +  (I-b)  T*^(l,b-I) 


+ 


-{ d-*2)  r . 


'+*'k(»-l,b)  +  (l-B)  cr*‘h 


*(«yJ‘hP(.-,,b)  -  oJ;hp(!-l,b)')  ]  } 

-  —  {  (l-X  )  [  B  (I  +1  ,b)  +  (l-B)  <yt"h(l+l,b)  1 

A»  '  J  L  *  X  J 


+  *3  [  8  cj*k(l+l,b)  +  (l-B)  <£h(l+l,b) 


(3.27) 
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The  stresses 


( i  •  I  ,b  ) ,  t1  (i-l,b).  <j*(i+l,b),  and  cr^fi-l.b),  which  are 
xy  xy  x  x 

required  in  the  equations  of7  motions  tor  the  mass  point  (i,b),  are  similarly 
determined  t rom  a  consideration  of  the  transmission  of  the  D'Alembert  forces 
on  mass  points  (i*l,b-l)  and  (i-l,b-l). 
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IV.  ILLUSTRATIVE  NUMERICAL  CALCULATIONS 


Numerical  calculations  for  a  large  number  of  test  problems  were 
performed  in  the  course  of  this  study;  these  were  necessary  to  verify  one  or 
more  aspects  of  the  capabilities  and  requirements  of  the  transmitting  boundary. 
Many  of  these  calculations  were  also  required  for  rejecting  or  modifying  a 
particular  conceptual  scheme.  In  particular,  all  the  numerical  schemes  described 
in  Chapter  III  were  the  results  of  extensive  conputat ional  experiments  and 
numerical  verifications.  Typical  among  the  problems  used  for  these  purposes 
are  those  described  in  the  sequel. 

A.  I  One  Dimensional  Propagation 

In  Sect.  3.3,  it  was  stated  that  the  proposed  transmitting  boundary  is 
theoretically  exact  for  linearly  clastic  propagation  in  one  dimension;  exact  in 
the  sense  that  no  additional  approximations,  other  than  those  associated  with  the 
basic  finite  difference  scheme,  arc  introduced  in  the  transmitting  boundary. 

This  is  verified  with  the  numerical  results  for  the  following  problems: 

(a)  A  semi -inf ini tc  one-dimensional  bar  subjected  to  a  sustained 
pressure  applied  at  the  fine  end.  The  exact  solution  for  this  case  is  known; 
i.e.,  the  applied  pulse  is  maintained  at  all  sections  along  the  bar.  The 
corresponding  numerical  solutions  obtained  with  Eq.  (3)  at  the  boundary  arc 
compared  with  the  exact  solution  in  Fig.  4.1.  Complete  time  histories  arc 
presented. 

(b)  A  second  problem  is  a  half-piano  subjected  to  the  pressure  pulse 
shown  in  Fig.  4.2,  applied  throughout  the  surface  of  the  half-plane.  In  this 
case,  the  pressure  pulse  consisted  of  both  loading  and  unloading  phases  of  the 
applied  pressures.  The  exact  solution  is  also  known,  which  is  the  applied  pulse 
reproduced  at  all  depths.  Numerical  results  calculated  using  the  transmitting 
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aundury  in  plane  strain  (i.e.,  Eqs.  3.7  through  3.15)  are  also  shown  in 
Fig.  4.2. 

The  r -suits  for  a  similar  proolem  in  which  the  pressure  pulse  is 
rep.aced  ay  a  shear  pulse  are  shown  in  Fig.  4.3.  The  calculations  were  also 
performed  with  the  plane  strain  model. 

(c)  Finally,  a  hal  f-space  is  subjected  to  a  normai  pressure  pulse 
app!  i  -d  uniformly  on  the  surface  of  the  half-space.  The  problem  was  formulated 
in  two  dimensions  under  ax i -symmetric  condition  fi.e.,  using  Eqs.  3.18  through 
3.2a);  however,  under  the  applied  uniform  surface  loading,  the  calculated  results 
should  be  the  same  as  those  for  plane  one-dimensional  propagation.  The  results 
shown  in  Fig.  4.4  verify  this. 

It  might  be  emphasized  from  these  results  that  there  are  no  reflections 
from  the  transmitting  boundary;  evidence  that  the  proposed  transmitting  boundary 
does  simulate  a  semi -I  if Ini te  (or  infinite)  space  "exactly"  in  the  case  of  plane 
propagation.  The  calculations  wore  all  performed  without  any  artificial  viscosity. 

■^•2  Plane-Strain  Cjlcu'ations 

Extensive  cal  cul  at  ional  experiments  were  also  performed  to  test  and 
verify  the  capabilities  of  the  transmitting  boundary  in  two  space  dimensions, 
under  both  piane  strain  and  axi-svmmetric  conditions. 

The  test  problems  in  plane  strain  generally  involve  a  partially  loaded 
half-plane,  as  shown  in  Fig.  415  or  4.9;  the  applied  loading  is  3  general  pressure 
pulse  described  in  Fig.  4.6  or  4.10,  whose  time  function  includes  both  loading 
and  unloading  phases  and  contains  a  wide  spectrum  of  frequencies.  The  problem  is 
i  sufficiently  general  test  of  the  capabilities  of  the  transmitting  boundary  in 
two  dimensions;  with  the  partial  loading  on  the  surface  of  the  half  plane,  truly 
two  dimensional  wave  effects  are  induced  at  the  boundary.  The  validity  of  the 
transmitting  boundary  is  verified  by  performing  two  sets  of  calculations 
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corresponding  to  two  different  depths  of  the  terminating  boundary.  Evaluation 
of  the  results  calculated  with  the  shallower  boundary  relative  to  those  obtained 
with  the  deeper  boundary  gives  a  means  of  verifying  the  capability  of  the 
transmitting  boundary. 

A  large  number  of  calcul ational  problems  were  performed;  problems 
involving  both  slow  and  fast  materials  were  considered.  However,  for  the  purpose 
of  illustrating  the  effectiveness  and  limitations  of  the  boundary,  the  results  of 
typical  cases  only  will  be  presented.  In  all  cases,  complete  time  histories  arc 
presented  for  stress  points  close  to  a  transmitting  boundary. 

The  half-plane  described  in  Fig.  4.5  has  a  slow  material  with  a  dilata- 
tional  velocity  of  1,600  fps.  The  transmitting  boundary  was  placed  at  two  depths; 
in  one  case  at  IS  ft.,  and  in  the  other  at  25  ft.  The  calculations  were  performed 
with  a  rectangular  grid  of  Ax  »  10  ft.  and  Ay  -  5  ft.,  and  a  uniform  time  increment 
of  At  *  0.75  msec,  was  used  In  the  numerical  integration.  This  represents  a  very 
severe  test  of  the  capability  of  the  transmitting  boundary;  observe  that  there  are 
only  3  mesh  lengths  to  the  boundary  in  one  case  (and  5  in  the  other),  such  that 
the  travel  time  from  the  surface  to  the  boundary  Is  only  about  1/6  the  duration 
of  the  applied  pulse  (and  about  1/4  the  duration  of  the  pulse  with  the  deeper 
boundary) . 

The  results  arc  presented  in  the  form  of  time  histories  of  vertical 
stresses  at  the  center  line  and  at  10  ft.  horizontally  from  the  center  line  (i.e., 
points  a  and  b  in  Fig.  4.5)  In  Figs.  4.7  and  4.8,  respectively.  It  should  be 
emphasized  that  the  stresses  presented  in  these  figures  arc  for  the  stress  points 
immediately  adjacent  to  the  transmitting  boundary  when  it  is  placed  at  the  shallow 
location  (i.e.,  0  *  15  ft.).  In  each  figure,  the  results  calculated  with  the 
transmitting  boundary  located  at  depths  of  15  ft.  and  25  ft.  arc  presented 
together  for  comparison;  on  this  basis,  the  validity  of  the  transmitting  boundary 
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is  clearly  evident  for  th:s  'aso. 

Plane  strain  ca'cul  at  ions  were  performed  also  for  a  partially  loaded 
half  plane  with  dilantional  velocities  of  3,100  fos  and  10,000  fps,  resembling 
soi  t  md  herd  rocks  r.-spect ’ vel y .  The  test  problem  is  idealized  as  shown  in 
fig.  1.9,  subjected  to  the  normal  pressure  pulse  shown  in  Fig.  4.13. 

Calculations  for  the  3,130  fps  mat:  rial  were  quite  we  1 1  -behaved ,  that  is, 
the  results  /ere  similar  to  those  for  the  1,600  fps  material  illustrated  earlier, 
and  consequently  will  not  be  repeated.  It  will  suffice  only  to  emphasize  that 
the  caltu'ations  were  performed  without  artificial  viscosity  and  serve  to  further 
verify  the  validity  of  the  transmitting  boundary. 

For  the  10,000  fps  material,  the  vertical  stresses  at  five  depths 
adjacent  to  the  center  line,  calculated  without  artificial  viscos'ty  (l.e.,  T  ■«  0)  , 
are  shown  in  Fig.  4.11.  Fig.  4.12  also  presents  the  vertical  stresses  at  roughly 
the  same  depths  and  at  a  horizontal  distance  of  120  ft.  from  the  center  line. 

The  results  presented  ii  these  figures  lllestrats  very  succinctly  the  difficulty 
that  could  arise  with  the  transmitting  boundary  when  a  fast  material  I*  involved. 
Although  the  main  pulses  are  correctly  transmitted  through  the  boundary,  high- 
frequency  oscillations  begin  to  appear  upon  un'oading;  these  oscillations  can  be 
quite  spurious  is  ill  as  t  rated  in  Figs.  4.11  and  4.12.  Such  oscillations  appear 
to  start  at  the  center  lime.  However,  the  fact  that  such  oscillations  are  of 
hi frequencies  suggest  that  artificial  viscosity  should  be  affective  in 
suppress :ng  then.  The  amount  of  viscosity  required  appears  to  depend  on  the 
propagation  velocity  of  the  material. 

in  Figs.  4.13  through  4.18  are  presented  the  calculational  results 
obtained  with  several  values  of  the  artlficia1  viscosity  coefficient  T,  combined 
with  two  different  time  increments.  From  the  results  given  In  Figs.  4.17  and 
4.18,  it  appears  that  for  a  material  with  10,000  fps  dilatational  velocity,  an 
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artificial  viscosity  of  T  =  0.40  is  sufficient  to  suppress  the  spurious 
oscillations  at  the  transmitting  boundary  without  affecting  the  main  pulses 
significantly.  Ail  the  calculations  were  performed  with  Ax  =  Ay  =  3b  ft. 

The  results  shown  in  Fig.  4.13  through  4.16  illustrate  the  consequence 
of  insufficient  viscosity  (r  =  0.28).  These  also  indicate  that  a  reduction  in 
the  integration  time  step  (from  1.30  msec,  to  1.40  msec.)  tends  to  reduce  the 
oscillations  but  does  not  eliminate  them. 

On  the  basis  of  the  two-dimensional  plane  strain  calculations  illus¬ 
trated  above,  the  following  observations  emc-rge: 

(i)  For  slow  material  (sa/  crf  <  3,000  fps),  the  transmitting  boundary 
can  be  used  for  plane  strain  problems  without  artificial  viscosity. 

(ii)  For  fast  material  (say  c^  >  3,000  fps) ,  the  transmitting  boundary 
ought  to  be  used  with  some  artificial  viscosity;  for  example,  with  c^  *=  10,000 
fps,  T  «  0.40  appears  to  be  appropriate  —  spurious  oscillations  are  suppressed 
without  affecting  the  real  signals  propagating  through  the  boundary. 

4 . 3  Axl -symmetric  Calculations 

The  transmitting  boundary  was  also  tested  for  a  variety  of  problems 
under  axially  f^mietrlc  conditions;  these  were  all  conducted  involving  an  axi- 
symmetric  half-space  subjected  to  normal  pressures  applied  at  the  surface, 
including  the  pulse  loading  shown  in  Fig.  4.13  and  a  periodic  sinusoidal  load 
shown  in  Fig.  4.26  applied  over  a  finite  circular  area. 

For  a  homogeneous  half-space  with  ■*  1,600  fps,  calculations  were 
performed  for  the  problem  shown  in  Fig.  4.20.  A  space  mesh  of  Ar  ■  10  ft.  and 
Az  *  5  ft.,  and  a  time  mesh  of  At  *  1.9  msec.,  were  used  In  the  calculations. 
Again,  in  order  to  provide  a  basis  for  judging  the  validity  of  the  results, 
the  calculations  were  conducted  with  the  transmitting  boundary  placed  at  two 
different  depths;  namely,  at  15  ft.  and  25  ft.,  respectively.  On  the  basis 


it  other  test  cas -s  .  tie  problem  presented  here  represents  the  most  severe  test 
ol  the  transit  1 1  i  ig  boundary  in  ax  i -syntiietr  ic  condition. 


The  vertical  stresses  at  a  depth  of  12.5  ft.  (stress  points  i mediately 
elare  the  transmitting  boundary  at  15-ft.  depth)  are  presented  is  Figs.  4.21 
and  1.22  tar  the  poi  its  a  and  b  in  Fig.  4.20  which  are,  respectively,  on  the  a.tis 
ot  s  m.ietry  and  at  r  ■  10  ft.  These  were  obtained  with  an  artificial  viscosity 
o’  t  «  0.15  js  required  for  cal  cu  I  at  ional  stability  of  the  axi -symmetr I c 
differencing  scheme.  In  each  or  Figs.  4.21  and  4.22,  the  results  calculated 
with  the  transmitting  boundary  located  at  15  ft.  and  25  ft.  are  presented  for 
the  pu-pose  of  comparison. 

A  similar  problem  involving  a  fast  materia',  crf  ■  10,000  f  ps ,  was  also 
analyzed.  The  problem  was  formulated  far  the  half-space  shoi/n  in  Fig.  4.23 
s  ib^ected  to  the  sale  pressure  pu'se  of  Fig.  4.19.  The  transmitting  boundary 
was  placed  at  a  depth  o*’  ISO  ft.,  and  calculations  were  performed  with  &r  «  Az  * 

•40  ft.  and  a  uniform  time  tesh  of  At  *  0.75  msec.  An  artificial  viscosity  of 
T  »  0.4  was  used. 

The  corresponding  results  for  the  time  histories  of  vertical  displace¬ 
ments  and  stresses  at  two  points  along  the  center  line  are  shown  In  Figs.  4.24 
and  4.25,  respectively.  It  might  be  emphasized  that  the  calculations  along  the 
center  line  were  found  to  be  most  critical;  i .e. ,  if  any  spurious  oscillations 
occur,  t  ley  in-ariab'y  start  from  the  center  line. 

The  transmitting  boundary  with  axial  symmetry  was  also  tested  for 
re peated  loadings.  For  this  purpose,  the  ha'  f-space  shown  In  Fig.  4.20  was 
subjected  to  a  sinusoidal  loading  function,  described  in  Fig.  4.26,  with  a 
period  of  120  msec.  The  vertical  stress  histories  are  presented  in  Figs.  4.27 
and  4.28,  respectively,  for  poi  its  a  and  b  of  Fig.  4.20.  In  each  figure,  the  calcu¬ 
lations  were  performed  with  the  transmitting  boundary  located  at  the  two  depths; 


namely,  at  IS  ft.  and  25  ft. 

It  is  of  significance  that  the  results  obtained  wi th  the  shal low  and 
deep  terminations,  as  shown  in  Figs.  4.21  and  4.22  and  in  Figs.  4.27  and  4.28, 
ara  very  close  to  each  other;  moreover,  the  calculations  remained  stable  long 
a'tcr  unloading  in  the  case  of  the  single  pulse  load,  whereas  in  the  case  of 
tie  periodic  load  the  sinusoids  are  accurately  reproduced.  From  the  experience 
with  several  other  schemes  tested,  it  is  the  general  observation  that  any 
transmitting  boundary  must  be  capable  of  virtually  perfect  simulation  of  toe 
effect  of  an  iifinitc  region,  otherwise  grossly  erroneous  calculations  occur 
almost  immediately  on  the  arrival  of  a  wave.  Therefore,  the  fact  that  the  main 
signals  arc  reproduced  correctly  through  the  transmitting  boundary  and  no  spurious 
oscillations  were  observed  on  unloading.  Is  strong  evidence  of  the  validity  of 
the  transmitting  boundary. 

4.4  Axl -symmetric  Layered  Elastic-Plastic  S/stem 

In  order  to  illustrate  the  capabilities  of  the  transmitting  boundary 
for  problems  Involving  complex  wave  forms,  calculations  were  performed  for  a 
3-laycr  elastic-plastic  axi-symmetr ic  half-space  shown  in  Fig.  4.29.  The  pressure 
pulse  applied  at  the  surface  of  the  half-space  is  the  same  as  that  shown  in  Fig. 
4.19.  The  three  layers  have  dilatational  speeds  as  follows: 

C|  ■  2,500  fps 
Cj  ■  6,000  fps 
Cj  *  10,000  fps 

The  calculations  were  performed  with  the  following  mesh  sizes: 

Ar  ■  40  ft.  (throughout) 


tz 

A  z 
A  z 


1 

2 
3 


20  ft. 
30  ft. 


10  ft 
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and  an  integration  time  slip  of  At  =  0.75  .nsec,  wa*;  used.  The  calculations 
were  done  with  artificial  viscosities  of  T  =  0.05,  0.15,  and  0.4  for  the  three 
respective  layers.  The  transmitting  boundary  was  placed  in  the  third  layer  at 
a  d'jp^h  of  180  ft. 

Two  sets  of  calculations  wore  performed.  In  one  case,  the  material  in 
all  the  layers  ./ere  assumed  to  be  elastic;  in  the  second  case,  the  material  In 
each  layer  was  assumed  ;o  be  el  asti :-perfectl y  plastic  with  a  yield  strength  of 
0.3  ks i ,  such  that  plastic  yielding  extends  all  the  way  to  the  transmitting 
boundary  as  depicted  in  Fig.  4.30. 

Typical  stress  histories  are  presented  in  Figs.  4.31  and  4.32  for 
points  on  the  axis  of  symmetry  and  at  depths  of  50  ft.  and  160  ft.,  respectively. 

In  both  of  these  figures,  the  elastic  and  elastic-plastic  stresses  are  presented 
for  comparison.  Similarly,  the  displacement  functions  for  the  elastic  and  elastic- 
plastic  calculations  are  shown  in  Figs.  4.33  and  4.34,  respectively,  at  depths  of 
50  and  160  ft. 

Again,  it  is  significant  that  the  results  are  quite  regular  with  no 
spurious  oscillations.  As  should  be  expected,  plastic  flow  tands  to  decrease  the 
maximum  stresses  at  a  stress  point,  but  larger  displacements  are  produced  under 
the  same  loading. 

From  this  last  problem,  two  points  should  be  emphasized: 

(I)  The  transmitting  boundary  is  capable  of  propagating  complex  waves 
such  as  those  arising  from  a  layered  system; 

(ii)  The  uffects  of  el astlc-plastlc  material  can  be  included  in  the 
transmitting  boundary. 

4 . 5  Results  of  lone  Sustained  Calculatioqal  Time 

To  illustrate  the  capability  of  the  transmitting  boundary  for  problems 
requiring  long  calcu’ atlonal  time  (i.o.,  many  calculation*!  time  slaps),  the 
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problem  described  in  Fig.  4.35  was  used.  It  resembles  crudely  a  problem  of  a 
nuclear  air  burst.  The  problem  is  Idealized  as  an  axi -symmetric  half-space 
subjected  to  an  initial  air  slap  applied  over  a  circular  region  of  radius  645  ft., 
plus  an  expanding  air  blast  whose  peak  pressures  decay  exponential!/  with  the 
radial  distance.  The  initial  air  slap  was  assumed  to  have  the  pulse  shown  in 
Fig.  4-10,  whereas  the  subsequent  air  pressures  correspond  to  those  described  in 
Ref.  [5]  for  a  one-megaton  surface  nuclear  burst.  The  half-space  consists  of 
three  layers  as  shown  In  Fig.  4.35  with  dilatational  velocities  of  2,500,  6,000, 
and  10,000  fps,  respectively.  The  thickness  of  layers  1  and  2  are  75  ft.  and 
160  ft.,  respectively;  whereas  the  third  layer  is  terminated  by  the  transmitting 
boundary  at  a  depth  of  345  ft.  from  the  surface.  The  material  in  each  la/er  is 
assumed  to  be  elastic-perfectl  y  plastic  with  a  uniform  yield  strength  of  10  ksi. 

The  calculations  were  performed  with  the  following  space  meshes: 

Ar  -  40  ft.,  for  all  layers 

Lz  -  10  ft.  for  layer  1 

Az  *  20  ft.  for  layer  2 

Lz  *  30  ft.  for  layer  3 

A  uniform  time  mesh  of  At  ■  I  msec,  was  used  throughout.  The  results  for  a 
duration  of  2.5  sec.  were  calculated;  corresponding,  therefore,  to  2,500  time 
steps,  representing  a  long  calculational  time. 

The  results  for  points  close  to  the  transmitting  boundary,  specifically 
at  a  depth  of  330  ft.,  ara  presented  in  Figs.  4.36  through  4.53.  These  arc 
complete  time  histories  of  motions  and  stresses  for  two  points  adjacent  to  the 
transmitting  boundary;  Figs.  4.36  through  4.44  are  at  a  radial  distance  of  1480  ft., 
vdicreas  Figs.  4.45  through  4.53  show  the  motions  and  stresses  at  a  radial  distance 
of  1320  ft. 
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These  results  clearly  sh ow  t'ae  capabilities  of  the  transmitting 
boundary;  specifically,  the  following  should  be  observed: 

(1)  There  are  no  spurious  oscillations,  even  after  long  calcul  ational 
time.  In  the  various  test  problems,  the  transmitting  boundary  has  been  shown  to 
transmit  correctly  the  main  pu'ses.  Moreover,  we  have  Indicated  that  the  boundary 
is  sensitive  to  small  inaccuracies,  such  that  unless  the  transmitting  boundary 
accurately  reproduces  the  effect  of  an  infinite  space,  the  errors  are  rapidly 
magnified  and  exhibited  in  the  form  of  spurious  oscillations. 

Therefore,  the  absence  of  spurious  oscillations  in  the  above  calculations 
is  an  indication  of  the  validity  of  the  results. 

(2)  The  calculated  particle  motions,  specifically  accelerations, 
velocitias  and  displacements,  are  essentially  zero  after  some  finite  time, 
reflecting  the  fact  that  the  particles  come  to  rest  after  the  passage  of  the 
stress  waves.  The  durations  of  the  signals  are  longer  for  points  at  larger 
distances  from  the  centar. 

(3)  The  stress  waves  are  quite  complex,  reflecting  the  effects  of  the 
material  layerings.  Aftir  the  passage  of  the  main  signals,  the  amplitudes  of  the 
stresses  are  significantly  reduced  (especially  at  the  1480*  range  as  shown  In 
Figs.  4.36  through  4  44).  This  Is  further  evidence  that  the  complex  waves  are 
being  transmitted  through  the  boundary  without  irregularities,  which  would  be 
expected  if  the  transmissions  through  the  boundary  were  inaccurate. 
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V.  OTHER  EXPLORATORY  STUDIES 


5. 1  Finite  Element  Type  Formulation 

The  notion  of  the  transmission  of  D'Alembert  forces  is  seemingly 
applicable  to  the  formulation  of  a  finite-element  type  transmitting  boundary. 

The  feasibility  of  such  a  formulation  was  explored.  For  this  purpose,  the 
plane-strain  model  shown  in  Fig.  5.1  was  used.  This  is  the  same  model  of  Fig. 

2.1  with  the  positions  of  the  elements  rearranged  as  shown. 

A  stress  point  and  the  four  neighborin  g  mass  points  may  be  considered 
to  constitute  a  finite  element;  this  would  be  equivalent  to  a  constant-strain 
rectangular  element.  A  stiffness  matrix  can  be  derived  for  such  an  element  and 
the  equations  of  motion  formulated  in  the  x  and  y  global  directions.  Calculations 
were  performed  on  this  basis  for  the  problem  described  in  Fig.  4.9,  in  which 
the  transmission  of  D'Alembert  forces  are  imposed  on  the  boundary  mass  points. 

In  Figs.  5.2  and  5.3  are  shown  the  results  for  the  vertical  and  shear 
stresses  at  several  depths  along  the  center  line.  False  signals  from  the 
boundary  appear  to  be  quite  evident  in  these  results.  These  error  signals  are 
slowly  oscillatory  and  grow  continuously.  In  contrast  to  the  formulation 
suggested  in  Chapter  III  in  which  artificial  viscosity  is  effective  in  suppressing 
the  high-frequency  error  growth,  there  is  no  way  that  the  type  of  errors 
illustrated  in  Figs.  5.2  and  5.3  can  be  easily  removed. 

On  the  basis  of  the  exploratory  calculations  performed,  as  described 
above,  the  proposed  transmitting  boundary  is  not  directly  applicable  to  the 
f Ini te-element  type  formulation.  At  least  some  modifications  of  the  schemes 
presented  in  Chapters  III  and  IV  will  be  required. 

5.2  An  Al  ternato  Approach 

Other  means  of  treating  the  effects  of  an  infinite  or  semi-inf  ini tj 
space  in  a  discrete-variable  solution  of  wave  propagations  nay  be  possible. 
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One  such  alternate  method  was  developed  as  described  In  the  Appendix.  Instead 
of  terminating  artificially  an  infinite  space  domain  with  a  transmitting  boundary, 
t'le  effects  of  the  infinite  space  may  be  included  in  the  discrete-variable 
formulation  within  a  finite  space  by  recognizing  certain  mathematical  properties 
of  wave  propagation  theory. 

The  basis  of  such  a  mathematical  scheme  Is  developed  in  the  Appendix 
for  a  centered  finite  difference  formulation.  Specific  finite  difference 
equations  (or  operators)  were  developed  for  orie-dimensiona!  and  two-dimensional 
plane  strain  el astic  medi a.  The  theoretical  basis  of  this  approach  is  described 
to  illustrate  an  alternate  di sere te-variabl a  treatment  of  an  Infinite  or  semi- 


V I .  SUMMARY  AND  CONCLUSIONS 


A  comprehensive  method  for  the  numerical  calculation  of  wave 
propagation  in  solid  media  is  summarized.  This  includes  a  physical  basis  for 
the  discretization  of  the  space  domain  that  leads  mathematically  to  a  centered 
finite  difference  system.  The  time  domain  Is  discretized  through  a  second-order 
step-wise  numerical  Integrator.  The  resulting  space-time  formulation  is 
equivalent  to  a  system  of  finite  difference  equations  for  a  set  of  discretized 
space  and  time  variables. 

It  might  be  emphasized  that  for  problems  involving  hyperbolic  or 
parabolic  systems,  such  as  in  dynamic  wave  propagation,  a  discrete-variable 
formulation  does  not  necessarily  imply  a  valid  numerical  solution.  Certain 
mathematical  conditions  for  stability  and  convergence  are  required,  which 
depend  on  the  resulting  set  of  equations.  In  the  formulation  suggested  herein, 
the  necessary  conditions  for  convergence  and  stability  are  available  in  the 
mathematical  theory  of  finite  differences.  On  this  basis,  explicit  stability 
conditions  for  plane  and  ax i -symmetric  propagations  were  developed. 

Calculations  with  any  discrete-variable  solution  method  are  necessarily 
limited  to  some  finite  space  domali.  For  problems  involving  Infinite  or  semi- 
infinite  spaces  the  "finite"  space  domain,  therefore,  must  be  terminated  in 
such  a  way  that  the  artificial  boundaries  of  the  pertinent  domain  will  reproduce 
the  effects  of  the  Infinite  space.  Such  a  boundary,  called  a  "transmitting 
boundary",  has  been  developed  for  the  cal cul at ional  method  described  above. 

The  proposed  transmitting  boundary  is  based  theoretically  on  the  concept  of 
the  step-wise  transmission  of  D'Alembert  forv.es.  it  simulates  exactly  the 
effects  of  an  infinite  space  if  the  transmission  speed  is  known.  Therefore, 
for  I i tear  one-di mens ional  clastic  propagation  (including  spherically  symmetric 
propagation),  the  proposed  transmitting  boundary  is  exact,  in  the  sense  that 
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no  errors,  other  then  those  associated  with  the  basic  finite  difference 
approximations,  are  introduced  aL  the  artificial  boundary.  This  concept  can 
be  readily  extended  to  propagation  problems  in  higher  dimensional  media,  or  to 
inelastic  ledia.  In  higher  dimensional  space,  the  formulation  of  the  necessary 
cal  cul  at ional  a'gorithm  at  the  boundary  is  facilitated  with  the  use  of  the 
appropriate  lumped-parameter  model;  the  step-wise  change  of  the  D'Alembert 
forc-s  on  a  mass  point  constitutes  the  underlying  basis  of  the  boundary.  In 
these  cases,  there  are  invariably  more  than  one  material  speed,  such  that  the 
actual  speed  of  transmission  of  the  D'Alembert  forces  is  not  known;  however, 
the  different  propagation  speeds  must  be  reconciled.  Among  the  several  schemes 
that  have  been  examined  for  this  purpose,  a  means  of  including  properly  all  the 
material  speeds  has  been  found.  The  resu’tlng  transmitting  boundary  was  tested 
thoroughly  with  a  large  variety  of  problems  and  cal cul ational  experiments.  Some 
approximation  in  the  transmitting  boundary  is  unavoidable  when  there  Is  more  than 
one  propagation  velocity.  The  errors  involved  in  this  approximation  invariably 
show  up  in  the  form  of  high-frequency  spurious  oscillations,  especia'ly  for  fast 
material,  which  normally  occur  after  passage  of  the  main  signals  (i.e.,  after 
unloading).  These  oscillations  are  characteristic  of  rapid  error  growth  which 
can  be  suppressed  through  .he  introduction  of  small  artificial  viscosity,  as  is 
normally  required  also  in  av.  I -symmetric  calculations,  which  does  not  significantly 
affect  the  real  waves. 

Based  on  the  studies  and  Investigations  conducted  herein,  the  following 
conclusions  may  be  emphasized: 

(I)  The  concept  of  the  step-wise  transmission  of  D'Alambcrt  forces  Is 
a  suitable  basis  fur  an  artificial  transmitting  boundary  to  simulate  the  effects 
of  an  infinite  or  seni-infinitc  space  in  one  and  two-space  dimensions. 
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(2)  For  problems  with  multiple  propagation  speeds,  it  is  of 


significance  to  observe  that  the  errors  associated  with  the  approximation 
required  to  reconcile  the  multiple  speeds,  have  a  tendency  to  grow  rapidly  into 
spurious  high-frequency  oscillations,  especially  for  fast  materials.  In  view 
of  this,  some  artificial  viscosity  should  be  used  with  the  transmitting  boundary. 
The  level  of  artificial  viscosity  normally  required  for  stability  of  axi -symmetric 
calculations  will  also  be  sufficient  to  suppress  the  error  growth  at  the  trans¬ 
mitting  boundary.  By  virtue  of  the  fact  that  the  boundary  is  extremely  sensitive 
to  small  errors  that  are  normally  exhibited  in  the  form  of  spurious  oscillations, 
the  absence  of  such  oscillations  in  a  particular  calculation  is  also  a  means  of 
ascertaining  or  verifying  the  validity  of  the  results. 

(3)  The  same  concept  is  applicable  for  a  transmitting  boundary  in 
three-space.  Although  no  difficulties  are  expected,  other  than  those  experienced 
and  resolved  In  the  two-space  problems,  the  generalization  of  the  basic  concept 
to  three-space  requires  additional  Investigation. 

(4)  Although  the  boundary  was  not  tested  for  any  formulation  in 
non-Cartesian  coordinates,  the  concept  of  step-wise  transmission  of  D'Alembert 
forces  is  applicable  to  any  coordinate  system.  Therefore,  no  special  difficulty 
is  anticipated  in  using  the  proposed  transmitting  boundary  for  teftoi.iat ing  a 
space  domain  with  an  irregular  geometry  based  on  the  generalized  did  .><*"*?• 
clement  model  described  in  Sect.  2.1. 

(3)  The  proposed  transmitting  boundary  is  also  applicable  to  complex 
wave  propagation  problems,  including  those  with  material  layerings,  periodic 
pulses,  and  other  realistic  loadings  such  as  blast  forces  and  earthquake 
disturbances. 

(6)  When  applied  to  clastic-plastic  material,  the  plastic  speed  of 
the  material  and  the  elastic  precursor  waves  must  be  reflected  in  the  formulation 
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of  the  transmi  ttin>j  boundary.  Although  the  extension  of  the  boundary  to 


elastic-plastic  Material  appears  ti>  be  valid, 
thoroughly  test'd  as  was  done  for  the  elastic 
concept  of  the  proposed  transmitting  boundary 
material,  if  the  proper  propagation  speeds  of 
of  doing  this  is  developed  herein. 


the  boundary  was  not  as 
cases.  Nevertheless,  the  basic 
is  suitable  for  elastic-plastic 
a  signal  are  included;  one  way 
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FIG.  2.21  AXISYMMETRIC  WAVE  PROPAGATION  -  RELATION  BETWEEN 
Ker  ANO  r  p*IO,  **0.00,  025,  0.40 


FIG.  2.22  AXISYMMETRIC  WAVE  PROPAGATION -RELATION  BETWEEN 
Ker  ANO  r  p»IS,  **0.00,  0.28, 0.40 
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FIG. 3.1  ONE  DIMENSIONAL  TRANSMITTING  BOUNDARY 


FIG.  3.2  TWO  DIMENSIONAL  TRANSMITTING  BOUNDARY 
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FIG.  3.3  TWO-DIMENSIONAL  MOOEL  FOR  PLANE  ELASTIC  MEDIA 
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FIG.  3.4  TRANSMITTING  BOUNOARY  FOR  AXISYMMETRlC  MEDIA 
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FIG.  4.7  VERTICAL  STRESSES  AT  X  «0,  Y-12.5  ♦» 
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FIG.  4.11  VERTICAL  STRESSES  REAR  CENTERLINE  (T«0) 


FIG.  4.12  VERTICAL  STRESSES  AT  X*I20  ft.  ( T  *0) 
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FIG.  4.16  VERTICAL  STRESSES  AT  X  -  120  ft.  ( I* «  0.28  ) 
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FIG.  4.17  VERTICAL  STRESSES  NEAR  CENTERLINE  ( I*  *  0.4  ) 
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FIG.  4.180  SHEAR  STRESSES  AT  X  »  I20»t.  CT  -0.40) 
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FIG. 4.33  VERTICAL  DISPLACEMENTS  AT  R*Oft,  i«90ft 
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FIG.  4.45  VERTICAL  ACCELERATION  AT  r>l920',  *030 
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FIG.  4.46  VERTICAL  VELOCITY  AT  r»l920\  *• 
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FIG.  4.48  RAOIAL  ACCELERATION  AT  r-1920',  Z -330 
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FIG.  4.49  RAOIAL  VELOCITY  AT  r»l920*, 
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RADIAL  DISPLACEMENT  AT  f!920\  I -MO* 
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FIG.  4.N  VERTICAL  STRESS  AT  r*l920’,  1  •  390 
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FIG.  4.52  RADIAL  STRESS  AT  r-1920*.  I  *330 
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FIG.  5. 
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FIG.  5.2  VERTICAL  STRESSES  AT  X  =120  f»  (  T  =0.28) 
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APPENDIX 


ALTERNATE  NUMERICAL  SIMULATION  OF  INFINITE  SPACE 


A.  I  General  Concepts 

The  proposed  alternate  numerical  procedure  Is  based  on  the  following 

the 


proposition:  Let  Un  ,  Wn  ,  Un  ,  and  Wn  be 

r  r  r,s;  u.  .  ’  r,s  u.  .  r,s;  u,  .  r,s  u.  . 

’  i  ,j  i  ij  '  I  ,j  i  ,j 

components  of  displacement  and  velocity,  respectively,  at  the  nth  time  step  and 

at  the  point  designated  by  r,s  due  to  a  unit  initial  displacement  at  point  i  ,j 

(u?  .  ■  I).  Similarly,  one  can  define  displacements  and  velocities  at  t  =  n  At 

and  the  point  r,s  due  to  u?  .  a  I ,  w?  .  *  I ,  and  w?  .  *  I. 

1  »J  •  »J  ■  *J 

For  linear  problems,  It  follows  that  the  components  of  displacement  and 
velocity  at  t  ■  n  At  and  at  the  point  r,s  are  given  by: 

“r.s  "  l.j  (uI,j  Ur,s;  uI(J  +  “l.j  Ur,s;  u,  }  +  WI,J  UI.s;  w,  }  + 


.  •  o  ,  ,n  \ 

+  WI.J  ur,ti*tiJ  1 


‘"r,.-  I,J  <  “°.j  . .  *?.«.  ♦ 


(A.  1 ) 


.  .•<>  i*.n  \ 

WI,J  Ur,s;  Wjj 


r.s 


f  .  (  u?  .  w”  +  i|  .  Iir  .  +  w?  ,  w"  t  + 

I  »J  I  ij  r.s,  I  »j  r,s,  U|  j  I  ,J  r,s,  Wj  j 


(A.2) 


,  .o  ,  ,n  \ 

WI»J  Ur,s;  Wj  j 


r,s 


1  (  u°  Wn  +  u°  Wn  .  +  w°  Wn  + 

l.j  {  i.J  Wr,s;  u,  j  I,j  Wr,s;  u,^  i,j  r,s;  w,  j 


(A.  3) 


+  w°  Wn  .  ) 

wl  ,j  Wr,s;  Wj  ^ 


(A.4) 
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where  the  surma t ion  extends  over  the  set  of  points  at  which  the  initial  conditions 
are  a I  I  non -zero. 

Therefore,  if  for  a  given  problem  the  quantities  Un  ,  , 

r,S*  Ui,j  r,S*  “i.j 

Wn  ,  Wn  ,  Un  . .  .  can  be  determined,  then  the  solution 

r,s;  u.  .  r,s;  u.  .  r,s;  u.  . 

'  .J  i  ,J  i  ,J 

of  the  discrete  ini t i al -val ue  problem  is  found  using  Eqs.  A.I  through  A.6.  For 

convenience  of  reference,  the  set  of  quantities  Un  ,  .  , 

r,s»  u|,j  r's*  u|,j 

Wn  ,  and  Wn  for  all  values  of  n  and  r,s  will  be  referred  to  as 

r,s;  u.  .  r,s;  u.  . 

'  *J  i  .J 

n  *p 

the  "operator  for  unit  initial  displacement  u  "•  Likewise,  U  .  ,  U  .  , 

r,S’  uf  ,j  r*‘*  ul,j 

Wn  and  Wn  .  for  all  n  and  r,s  constitute  the  "operator  for  unit 

r.s;  u,  j  r.s;  u|(J 

initial  velocity  u  In  the  more  general  case  discussed  above  which  actually 
corresponds  to  axi-symmetrlc  wave  motion,  one  has  to  compute  four  unit  operators; 
l.e.,  for  u,  u,  w,  and  w,  respectively. 

The  construction  of  the  unit  operators  is  the  stage  at  which  the  aetjal 
integration  of  the  dl f ferentia! -di f f-arence  equation  by  means  of  the  appropriate 
numerical  Integrator  is  performed.  The  construction  of  these  unit  operators  may 
present  some  difficulties  under  certain  circumstances;  for  certain  classes  of 
problems,  however,  some  simplifications  are  possible: 

I)  Due  to  the  numerical  integrator  used  (i.e.,  the  Newnark  8-method), 
a  unit  disturbance  at  a  point  for  t  ■  0  with  all  other  initial  conditions  being 
equal  to  zero  will  create  non-zero  displacements  or  velocities  only  at  a  finite 
number  of  points  at  t  *  &t;  the  unit  operator  will  propagate  In  a  pyramid-like 
fashion  into  the  time  dor, tain  with  the  apex  of  the  pyramid  located  at  t  ■  0. 

This  permits  considerable  savings  in  the  computations  necessary  to  construct  the 
finite  difference  operators. 
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2)  For  one-dimensional  and  plane  strain  wave  propagation  problems, 
if  central  finite-difference  approximations  are  used,  then  the  operators  will  be 
either  symmetrical  or  anti-symmetrical  with  respect  to  the  coordinate  axes  passing 
through  the  apex  of  the  pyramid.  The  di f ferential -di fference  equations  obtained 
on  the  basis  of  the  appropriate  lumped-parameter  models  are  central  finite- 
difference  analogues  of  the  corresponding  continuous  formulations.  In  this  case, 
a  substantial  economy  of  computations  is  realized. 

A. 2  One-Olmenslonal  Wave  Motion  -  Infinite  Spatial  Domain 

The  equation  of  motion  for  one-dimensional  wave  propagation  when 
formulated  using  the  lumped-parameter  model  of  Fig.  2. 1» 

U,  -  (—^-)2  (u,+2  -  2  Uj  +  U| _2>  .  (A. 5) 

Ax 

where  Ax  is  the  uniform  space  increment  in  the  x-directlon,  and  c^  is  the 
dilatatlonal  velocity  of  propagation.  The  domain  of  integration  is  assumed  to 
be  restricted  to  (see  Fig.  A.l) 


0  <  t  <  N  At 
0  <  x  <  p  Ax 


Both  at  x  ■  0  and  x  ■  p  Ax,  no  boundary  conditions  are  prescribed. 


The  unit  operators  to  be  constructed  are  those  for  u 


°  i  j  *o  , 

I  ■  I  and  u  j  ■  I . 


The  construction  is  made  by  integrating  Eq.  A. 5  using  the  quadrature  relations 
of  Eqs.  (19)  and  (20);  this  process  is  simplified  upon  realizing  the  symmetry 
of  the  operators  with  respect  to  the  apex  point  i.  in  other  words, 

Un  -  un 
l+r;  u,  Ui-r;  u, 


r,n  j*.n 

Ul+r;  u,  "  Ul-r;  u, 
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u?  . 

i  -4 ;  u . 


u . 

i 


(A.S) 


Using  this  symmetry,  the  expression  for  acceleration  at  the  axis  of  symmetry  of 
t'ie  unit  operators  becomes: 


=  2 


(ui+2 


(A.7) 


The  final  values  of  the  displacements  and  velocities  are  obtained  from 


Ip  /  o  ,,n  .  .o  „n  » 

1.1  (  ui  Ur;  u,  +  ul  ur;  u,  > 


f.,  «  “?  K:  u,  *  «?  K;  6,  > 


(A.8) 


In  order  to  realize  the  most  economical  use  of  the  procedure  outlined 
above,  three  distinct  situations  should  be  considered  (see  Fig.  A.I): 

(1)  N  <  p-l.  In  this  case,  the  operators  should  be  constructed  up  to 
and  including  the  time  level  n  ■  N.  The  triangular  construction  with  apex 
centered  at  the  extreme  ends  of  the  domain  of  Integration  will  suffice  since  it 
corresponds  to  the  most  critical  situation*  It  should  be  observed  that  the  points 
to  the  right  of  the  dotted  line  in  Fig.  A.I  are  unaffected  by  the  unit  effects  at 
the  apex  since  such  points  are  not  Inside  the  "influence  areal1  of  the  apex. 
Consequently,  the  integration  is  carried  out  for  considerably  less  points  than 
the  total  of  points  within  the  domain  of  integration. 

(2)  p-l  <  N  <  2 ( p- 1 ) .  The  pyramid-like  operators  already  computed  up 
to  and  Including  n  -  p-l  can  be  used  to  evaluate  the  values  of  displacements  and 
velocities  for  points  inside  the  spatial  domain  of  integration  for  all  time  stsps 

up  to  and  including  n  **  2 ( p— 1 ) •  Eq.  A.8  must  be  used  for  this  purpose. 
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In  extending  the  operators  up  to  the  time  step  n  ■  2(p-l),  no  data  pertaining  to 
points  outside  the  prescribed  domain  of  integration  are  computed. 

(3)  N  >2(p-l).  In  order  to  construct  the  unit  operators  to  a  depth 
of  n  =  N  >  2 ( p— I ) ,  the  triangular  portions  evaluated  in  (I)  above  for  up  to 
n  »  p-l  can  be  used  with  the  points  at  n  *  2(p-l).  However,  the  influence  area 
of  the  point  I  at  n  •  2 ( p- I )  goes  well  beyond  the  spatial  domain  by  an  amount 
equal  to  (p-l)  Ax.  Therefore,  In  order  to  carry  out  the  construction  beyond  n  » 
2(p-l),  p-l  additional  points  outside  the  domain  of  Integration  must  be  evaluated 
at  n  -  2 (p-l ) . 

The  number  of  extra  points  to  be  evaluated  In  order  to  carry  out  the 
integration  up  to  n4J  depends  upon  the  number  of  mats  points  Inside  the  domain  of 
integration.  The  obvious  advantage  of  the  algorithm  is  that  for  N  2  2 (p-l )  the 
total  number  of  points  which  are  integrated  are  less  than  the  sum  of  all  points 
inside  the  domain  of  integration;  and  beyond  2 (p-l )  additional  axtre  points  must 
be  integrated  for  only  at  intervals  of  p-l.  It  should  be  emphasized  that  after 
the  construction  of  the  operators,  the  extra  points  need  not  be  retained  in  the 
memory  of  the  computer.  Furthermore,  the  number  of  additional  points  to  be 
considered  at  each  time  step  of  multiples  of  p-l  after  n  ■  2(p-l)  increases  by 
an  amount  of  p-l  up  to  the  middle  of  the  Interval  [p-l,  N]  ;  it  then  decreases  by 
an  amount  of  p-l  gradually,  thereby  avoiding  superfluous  computations. 

Once  these  unit  operators  are  obtained,  it  is  a  straightforward  matter 
to  compute  the  displacements  and  velocities  for  a  prescribed  set  of  initial 
conditions  using  Eq.  A.8.  An  advantage  of  the  proposed  method  is  that  the 
operators  derived  for  a  given  grid  and  a  time  increment  can  be  used  to  solve 
different  wave  propagation  problems  corresponding  to  different  Initial  pulses. 
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A.  3  Plano  Strain  Wave  Motion  -  Infinite  Spatial  Domain 


The  discrete  equations  of  motion  for  plane  strain  wave  propagation 
formulat  'd  on  the  basis  of  the  plane-strain  I  jmped-parameter  model  are: 

t  J  J  C  M 

u.  .  =  ("T — )  (  u.  ,  .  -  2  u.  .  +  u..«  .  )  +  (—/■*)  (  u.  .  ,  - 

i  ,j  Ax  i  -2  ,j  i  ,j  i+2 ,  t  'Ay  '  '  I  ,j-2 


l+l.j+l  “  Wl+I,j-I 


-  w,  ,  ,  + 

i-l  ,j+l 


w. 


i-l.j-l 


(A.9) 


*!.j  ’  (-“’2  (  "l.J-2  '  2  "l.J  +  1  + 


*  2  “l.j  +  "|M,1  >  + 


c2-c2 

s,  Sy  +  (  “W.J+1  ‘  “'♦'J-l 


*  “l-I.J+l  +  Vl.j-I  > 

The  domain  of  integration  is  defined  by 


0  <  t  <  N  At 
0  £x  <  1/2  p  Ax 
0  <  y  <  1/2  q  Ay 


(A. 10) 


(A.II) 


The  spatial  domain  of  integration  is  actually  infinite;  Eq.  A.II  defines  a  portion 
of  it  over  which  numerical  integration  will  be  performed. 

It  is  observed  that  the  first  and  the  second  terms  in  Eqs.  A.9  and  A.  10 
are  symmetrical  and  the  third  terns  in  both  of  these  equations  are  anti-symmetrical 
with  respect  to  the  reference  coordinates  passing  through  the  point  I  ,J .  This 
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observation  leads  to  considerable  simplifications  in  the  construction  of  the 
unit  operators,  because: 

(1)  In  the  operators  for  unit  initial  displacement  u?  .  and  unit 

*  »J 

initial  velocity  u?  .  ,  Un  ,  Un  ,  u”  .  and  u”  .  are 

r*s;ul,j  r*s;ui,j  r*S;ui,j  r*s:ul,j 


,*,n 


symmetrical  and  W‘  ,  W  ,  W  and  W 

r's:ui,j  r*S:ui,j  r’s;ui,j  r'S;ui,j 

symmetrical.  For  example. 


are  anti* 


yO  ^  yll  ^  yfl  ^  yl) 

l+r,j+s;  u,  j  “  ui-r,j+s;  u.  .  "  l+r,j-s;  u,^  “  l-r,j-s;  u,^. 


Wn  m  uO  m  nH  B 

l+r,j+s;  u,tj  i-r,j+s;  u,^  l+r,j-s;  u^j  l-4,j-s;  uf  j 

(2)  In  the  operators  for  unit  initial  displacement  w?  ,,  and  unit 

1  iJ 


initial  velocity  w°  ,  uj  ,  if  uj  .  and  Oj  . 

ifJ  r.s,  w,tj  r.s,  w,(j,  r,s,  w|(j  r’*»w|,j 

ant  I -symmetrical  and  w"  .  w"  .  ,  w"  .  and  ll"  .  ar 

r,*'  ”,  j  r»*’  w|,j  r**’wl,j 


are 


symmetrical  with  respect  to  the  coordinate  axes  through  i,J. 

Therefore,  in  constructing  the  operators  for  u?  .  ■  I ,  u ?  .  ■  I , 

■  >J  ' 

w?  .  ■  I,  and  w,  ,  ■  I,  it  is  sufficient  to  integrate  only  for  points  which  are 
•  *J  ■  ,J 

included  inside  a  quadrant  of  the  reference  coordinates  originating  from  I  ,J 
(see  Fig.  A.2).  For  points  in  which  either  I  or  j  is  I  or  2,  the  computation  of 
accelerations  cannot  be  made  through  Eqs.  A. 9  or  A.  10  directly.  The  use  of 
symmetry  or  anti-svmmetry  yields  the  following  simpler  expressions  at  such  points: 


For  u«(,  -  I  or  u«f|  -  I 


1,1 


2  A*  (u3  ,  -  u  ,)  +  2  (-^-)2  (U)  ,  -  u‘>  + 
Ax  *'  Ay  '*3 
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(4.12a) 


2  2 
c  ,  •  c 

+  4  _ _ s_ 


Ax  Ay 


2,2 


'2,2  (  ,  'w2 ,4  "  3  w2,2*  +  ^  ^  (w4,2  "  3  W2,2*  + 

Ay  Ax 


c  o 

/  s  >2 


2  2 
Cd  ”  Cs 

+  - -  (u3  3  ""3  |  "u|  3  +  ul  |> 

Ax  Ay  3(5  3,1  '*3 


(A. 12b 


"l.i  '  2  (f’2  (u3.i  ■  Ui.i>  +  <^-'2  <“l.j-2  '  2“'.j  *  “l.j+2>  + 
2  _  1 

+  2  TT^  (W2*j+l  *  w2,j -i ) 

Ax  Ay 


(A. 12c) 


ui,i  *  <“^“)2  (ui-2,i  '2  “i.i  *2  (-f->2  *“l,»  ‘  “l.l*  + 


c2  -  c2 

*'1^7  <"'+l.2'”i-l,2> 


(A. 1 2d) 


:2.i  '  (-^)2  (”2.i-2  *  2  "2.J  +  l'2,J«>  +  <-^>2  <"4,J  ‘  3  »2,J>  + 


2  2 
c .  -  c 

+  - - (u,,  . , ,  -  u,  .  ,  -  u,  ...  +  u,  ,  |) 

Ax  Ay  3«J+I  3'J"1  ,»J+I  ’'J-1 


(A. 1 2e) 


!i,2  •*  (wt  ,4  *  3  w|  ,2}  +  (_^*)Z  (w!-2,2  ‘  2  wl,2  +  wi+2,2)  + 


2  2 
c.  -  c 

+  k  A  (ul+l,3  “  u?+l,1  "  ui-l,1  +  UM,1* 

Ax  Ay 


*1.1  *“2,2  "”l,j  "”1.1  -“2.J  **1,2  *° 


(A.  1 2  f ) 
(A.l  2g) 


And  for  w?  .  *  I  or  w?  .  *  I  : 
1,1  1.J 
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(A. 1 3a) 


«i,|  *  2  (wl,3  "  "l,l)  +  2  ("^“)2  (W3,I  *  “».l)  + 


2  2 
c.  "  c« 

+  4  — - -  u 

Ax  Ay 


2,2 


S  v2 


,2  *  (“}  (U4.2  *  3  U2,2}  +  <  ^  >  (U2,4  '  3  u2,2}  + 


2  -  c2 

+  7TT  '”3.3- “3.!  •«l.3t»M> 


“>,j  '  <_^*)2  <"l,J-2  '  2  "l,J  +  “l  .Jt-2*  +  2  (“3.J  *  "l.j*  + 


+  2^  (^.J+'*“2-J-') 


W,  ,  ■  2  (— ^-)2  (wl  3  -  w,  ,)  +  (— *-)  (w,.2  ,  •  2  W|.|  +  W1+2.P 


'1,1 


Ay 


2  2 


Ax 


i-2,1  ‘  "1,1  "1+2, r 


*2TTi  (Vl,2 

Ax  Ay 


u2,j  “  (^”)2  (u4,j  “  3  U2,J}  +  (u2 ,J ”2  ’  2  U2,J  +  U2,J+2 


c2  c2 

+ir^f  <W3.j+»  ■  w3.j-'  "wu+* +wu-i) 


>  + 


S  n2 


u|  ,2  *  <“>  (ui-2,2  ’  2  Ul,2  +  ul+2,2)  4  (~}  (ul  ,4  "  3  4 

(wl+J,3  "  wi-l ,3  '  wl+l,l  +  wl-l,l) 


2  2 
c .  -  c_ 


Ax  Ay 


(A. 13b) 


(A. 13c) 


(A. 13d) 


(A. 13e) 


(A. 1 3f) 
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-  0 


(A. 1 3  g) 


M  •• 


Furthermore,  it  should  bo  noted  that  the  operators  for  u®  |  ■  I  or 

ik  .  =  I,  ii.  .  =0  when  i  end  j  arc  both  even,  and  w.  .  ■  0  when  I  and  j  are 
■  •  •  •  »J  1  ij 

both  odd.  Conversely,  the  operators  for  w?  ,  ■  I  or  w?  ,  ■  I,  u,  .  “0  when  I 

1,1  i,i  •  »J 

end  j  are  both  odd,  and  w.  .  »  0  when  i  and  j  are  both  even.  This  is  a  result 

•  ,J 

of  the  forms  of  Eqs.  A. 9  and  A. 10. 

When  these  considerations  are  used,  the  construction  of  the  operators 
for  unit  displacements  and  velocities  becomes  appreciably  simpler. 

Depend  in j  upon  the  size  of  the  spatial  domain  of  integration  and  the 
number  of  time  steps  at  which  integration  is  desired,  three  different  situations 
are  possible  (see  Fig  A. 3).  If  it  is  arbitrarily  assuned  that  p  £  q,  then: 

(1)  N  <  1/2  q  -  I.  In  this  case,  the  number  of  points  at  vAich 
integration  is  to  be  performed  in  constructing  the  unit  operators  will  be  much 
less  than  the  total  number  of  points  inside  the  domain  of  integration.  The 
operators  for  unit  displacements  and  unit  velocities  extend  into  the  time  domain 
in  a  pyramid-like  fashion,  and  the  integration  must  be  carried  out  for  n  •  N  time 
steps. 

(2)  1/2  q  -  I  <  N  £  2  (1/2  q  -  1).  The  most  economical  way  of 
constructing  the  operators  Is  to  first  construct  the  pyramid-like  portions  p  to 
and  including  n  -  1/2  q  -  I,  and  then  to  use  these  operators  in  association  with 
Eqs.  A.  1  through  A. 4  for  the  computed  values  of  displacements  and  vclocitiss  at 
n  =  1.2  q  -  I;  in  this  manner  for  time  steps  beyond  1/2  q  -  I,  only  the  points 
inside  the  spatial  domain  can  be  computed  without  the  necessity  of  going  beyond 
the  imaginary  boundaries. 

(3)  N  >  2(1/2  q  -  i).  For  n  <  2(1/2  q  -  I),  the  construction  of  tha 
unit  operators  is  done  by  the  procedure  described  above.  However,  at  n  ■ 
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2(1/2  q  -  I)  the  wave  front  is  already  outside  the  spatial  domain  of  Integration; 
consequently.  In  order  to  compute  for  the  points  at  n  ■  2  (1/2  q  -  I),  some 
additional  points  outside  the  domain  of  integration  at  n  -  2(1/2  q  -  I)  must  be 
evaluated  (see  Fig.  A. 3).  The  information  about  the  displacements  and  velocities 
at  the  points  inside  the  triangular  area  ABC  Is  sufficient  to  construct  the 
operators,  by  the  aid  of  the  pyr»ld-llke  portions  of  the  operators,  up  to  n  ■ 
3(1/2  q  -  I)  without  having  to  compute  for  any  points  outside  the  spatial  domain 
of  integration.  In  this  way,  only  at  Intervals  of  1/2  q  -  I  beyond  n  -  2(1/2  q  -  I) 
is  it  necessary  to  consider  additional  points  outside  the  spatial  domain.  Again, 
following  the  construction  of  the  operators,  the  Information  concerning  the 
additional  points  outside  the  domain  of  Integration  need  not  be  retained  in  the 
memory  of  the  computer.  Finally,  the  number  of  extra  points  for  which  computations 
are  to  be  made  increases  gradually  In  the  Interval  (1/2  q  -  I,  N)  up  to  the  middle 
of  this  Interval,  and  then  may  be  decreased  gradually  to  only  the  required  points 
inside  the  domain  of  Integration  at  n  -  N. 

The  next  step  following  the  construction  of  the  unit  operators  is  the 
evaluation  of  the  displacements  and  velocities  of  the  mass  points  using  the 
principle  of  superposition  as  expressed  by  Eqs.  A.I  through  A.4.  Obviously,  the 
construction  of  the  unit  operators  depends  on  the  physical  properties  of  the 
material  in  the  domain  of  Integration,  the  mesh  size  and  the  time  increment,  but 
independent  of  the  Initial  conditions.  Consequently,  the  unit  operators  derived 
for  a  given  problem  can  be  used  for  different  sets  of  Initial  conditions. 

A. 4  Treatment  of  Semi-Inf  ini  te  Domain. 

The  algorithm  described  up  till  this  point  is  capable  of  handling 
infinite  spatial  domains  of  Integration  with  no  real  boundaries.  When  boundary 
conditions  are  also  prescribed  In  addition  to  the  initial  conditions,  the  proposed 
algorithm  can  easily  be  adapted  to  take  care  of  semi-infinite  domains.  The  modified 
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algorithm  presented  below  .s  for  one-dimensional  wave  propagation  In  a  seml- 
iiflnite  domain;  extension  to  plane  strain  wave  propagation  In  a  semi-inf  Ini te 
domain  is  s t ra i gh tforward . 

(1)  Let  AB  represent  the  real  boundary  In  the  (x,  t)-plane  where 
boundary  conditions  are  known  (see  Fig.  A. 4).  The  spatial  domain  of  Integration 
;s  limited  by  an  imaginary  boundary  at  the  right-hand  side  of  the  domain.  Let  us 

arbitrarily  suppose  that  the  Initial  conditions  are  given  at  points  2,  3,  . 

p',  where  p1  <  pH,  and  beyond  the  point  p*  Initial  conditions  are  uniformly  zero. 

(2)  The  operators  for  unit  displacement  and  unit  velocity  are 
constructed  for  n  =  p-2  In  the  same  way  as  that  described  for  an  Infinite  medium. 
Since  point  I  i>  the  boundary  point,  the  depth  of  the  operators  must  be  p-2  time 
steps,  and  not  p-l ,  which  was  the  case  previously. 

(3)  The  points  at  and  to  the  right  of  the  line  CD  up  to  the  time  step 
n  =  p-2  are  not  influenced  by  the  boundary  conditions  along  AB.  Consequently, 

the  displacements  and  velocities  for  such  points  can  be  computed  using  the  derived 
operators  and  Eq.  A. 8. 

(4)  At  each  time  step  between  n  *  I  and  n  ■  p-2,  the  points  to  the  left 
of  the  line  CD  are  affected  by  the  boundary  conditions  along  AB  as  well  as  the 
initial  conditions  at  t  =  0.  Therefore,  the  displacements  and  velocities  at  such 
points  cannot  be  computed  by  the  unit  operators.  However,  it  should  be  observed 
that  at  each  time  step  between  n  =  I  and  n  =  p-2,  the  portions  of  the  spatial 
domain  enclosed  by  the  lines  AB  and  CD  are  bounded  such  that  at  both  boundaries 
the  displacements  and  velocities  are  known.  Consequently,  the  points  between  AB 
and  CD  can  be  integrated  directly,  using  the  numerical  Integration  procedure  for 

a  bounded  domain. 

(5)  If  Integration  is  desired  for  n>  p-2,  then  the  above  steps  should 
be  repeated  for  p-2  <  n  <  2(p-2).  Again,  the  boundary  conditions  along  AB  do  not 
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have  any  effect  on  the  points  at  and  to  the  right  of  the  line  CD1;  such 
points  are  integrated  using  the  initial  conditions  at  n  *  p-2  with  the  unit 
operators.  The  points  to  the  left  of  CD'  are  Integrated  directly  as  before. 
However,  since  the  wave  front  shown  by  the  line  EFGH  is  beyond  the  spatial 
domain  at  n  =  p-2,  in  computing  for  points  for  n  £  p-2,  p •  -2  additional  points 
for  which  the  displacements  and  velocities  are  determined  while  computing  for 
0  <  n  <  p-2  must  be  taken  into  consideration.  Similarly,  if  N  >  2(p-2),  then 
at  n  *  2 (p-2) ,  an  additional  (p'-2)+(p-2)  points  must  be  computed. 

As  before,  the  number  of  additional  points  beyond  the  spatial  domain 
of  integration  that  are  to  be  computed  at  time  steps  which  arc  multiples  of  (p-2) 
increases  in  number  by  p-2  at  each  p-2  time  steps  up  to  the  middle  of  the  interval 
(p-2,  N)  and  then  gradually  decreases  to  the  size  of  the  spatial  domain  at  n  ■  N. 
This  is  the  most  economical  way  of  attacking  the  problem.  Also,  it  should  be 
realized  that,  while  computing  for  n  >  k(p-2),  the  extra  points  computed  at 
n  <  (k- 1 ) (p-2)  need  not  be  retained  inside  the  computer. 

The  algorithm  described  In  this  section  indicates  that  the  use  of  the 
proposed  numerical  procedure  is  not  limited  by  the  assumption  that  the  Initial 
impulse  must  be  restricted  to  within  the  spatial  domain  of  Integration;  as  long 
as  the  wave  front  Is  faithfully  followed,  the  disturbances  outside  the  limited 
spatial  domain  can  also  be  accounted  for. 

A. 5  Concluding  Remarks 

The  proposed  algorithm  Is  based  on  the  assumption  of  a  linear  problem 
and  the  applicability  of  the  principle  of  superposition.  It,  therefore,  cannot 
be  used  when  nonl  Inearl  ties  due  to  material  behavior  or  large  displacements  are 
present.  Also,  It  Is  limited  to  Newmark's  integrator  with  p  ■  0;  non-zero  values 
of  B  do  not  lend  themselves  easily  to  the  proposed  numerical  procedure  since  the 
unit  operators  will  not  be  as  simple  as  those  for  6  ■  0.  Other  numerical 
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integrators  would  also  have  the  same  complication. 

Despite  taese  limitations,  the  algorithm  is  promising  in  elastic  wave 
propagation  problems  in  one-dimensional  and  plane  strain  media.  It  can  readily 
be  generalized  to  include  wave  motion  In  spherically  symmetric  and  axl-symmetric 
media,  as  well,  since  the  only  change  comes  in  the  construction  of  the  unit 
operators  for  displacements  and  velocities. 

The  proposed  numerical  algorithm  has  the  following  advantages: 

(1)  In  addition  to  those  associated  with  the  finite-difference  scheme 
itself,  no  additional  approximations  are  introduced. 

(2)  The  numerical  procedure  Is  equivalent  to  the  solution  of  the 
finite-difference  scheme;  it  then  follows  that  the  stability  and  the  convergence 
properties  of  the  original  difference  scheme  remains  unaltered  by  the  proposed 
algorithm,  and  that  no  other  types  of  Instabilities  are  introduced  Into  the  problem. 

(3)  The  operators  for  unit  displacements  and  velocities  depend  on  the 
material  propertias  of  the  medium,  the  mesh  size  and  the  uniform  time  step. 
Therefore,  once  they  are  constructed,  they  can  be  used  to  solve  the  same  problem 
with  dlffarent  sets  of  initial  conditions. 

The  mala  drawback  of  the  procedure  stems  from  the  fact  that  depending 
on  the  number  of  discrete  mass  points  Included  in  the  spatial  domain  and  the 
number  of  time  steps  for  which  Integration  Is  desired,  It  may  be  necessary  to 
compute  for  some  additional  points  outside  the  spatial  domain  in  order  to  follow 
the  wave  front.  Ordinarily,  the  wave  front  must  be  Inevitably  followed  by 
computing  for  a  gradually  increasing  number  of  points  outside  the  spatial  domain 
of  integration  at  each  time  step;  the  proposed  algorithm  eliminates  this  by 
offering  a  compromise  in  that  It  is  necessary  to  consider  extra  points  only  at 
Intervals  of  the  time  increment  which  is  the  depth  of  the  unit  operators.  Thus, 
depending  on  thu  mesh  size  and  the  total  number  of  time  steps  N  for  which 
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integration  is  to  be  performed,  one  may  totally  avoid  the  extra  computations 
or  limit  the  additional  points  to  such  a  number  that  the  whole  solution  process 
becomes  no  less  economical  than  the  conventional  integration  procedure  which  is 
incapable  of  handling  the  problems  arising  at  the  boundaries.  For  any  given 
problem,  as  soon  as  the  problem  parameters  are  fixed,  the  analyst  may  easily 
determine  the  most  economical  mesh  configuration  for  which  a  minimum  of 
additional  points  is  necessary.  Table  I  Is  designed  to  illustrate  these 
concepts  in  the  case  of  wave  motion  In  a  one-dlmenslonal  medium  extending 
infinitely  in  both  directions. 

The  algorithm  has  been  successfully  applied  to  a  number  of  problems 
in  one-dimensional  and  plane  strain  media.  It  has  been  noted  that  the  construc¬ 
tion  of  the  unit  operators  can  be  done  with  a  relatively  high  efficiency  In  a 
digital  computer.  The  construction  of  these  operators  constitutes  the  bulk  of 
the  solution  process;  however,  for  a  finite  number  of  time  steps  which  Is  not 
excessively  large  compared  to  the  minimum  number  V  Aais  points  along  any  one 
of  the  spatial  directions,  this  can  be  economically  accomplished,  and  needs  to 
be  done  only  once  for  a  class  of  problems. 
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FIG.  A.  I  CONSTRUCTION  OF  UNIT  OPERATORS  I  ONE  -DIMENSIONAL 
PROPAGATION  IN  AN  INFINITE  SPACE 
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FIG.  A. 2  CONSTRUCTION  OF  UNIT  OPERATORS  I  PLANE  STRAIN 
PROPAGATION  IN  AN  INFINITE  SPACE 
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FIG.  A.4  INTEGRATION  FOR  ONE -DIMENSIONAL  WAVE  MOTION  :  IN  A  SEMI-INFINITE  SPACE 
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